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Foreword 


It is with pleasure that I write the foreword to this excellent book entitled “Fractal 
Solutions for Understanding Complex Systems in Earth Sciences” edited by 
Professor V.P. Dimri. The book deals with a role of fractals and multi-fractals in 
understanding problems of earth sciences and their solutions. The chapters on 
potential fields suggest that mono-fractals or homogeneous functions are scaling 
function, whereas multi-fractals and multi-homogeneous models have to be dealt 
separately by specific techniques. Some of the chapters deal with Detrended 
Fluctuation Analysis (DFA) for better understanding the intrinsic self-similarities 
for non-stationary well-log geophysical data. The Fractal Differential Adjacency 
Segregation (F-DAS) is found better than the conventional box counting technique 
to estimate fractal dimension. Amongst different types of earthquakes analysed, the 
multi-fractal properties are more pronounced for signals with distinct P, S and coda 
waves. Finally, an overview of fractal approach to fire problems is worth 
mentioning. 

The editor, Professor Dimri, is to be congratulated on putting together various 
aspects of fractal and multi-fractal in the problems of earth sciences and their 
solutions. The book is recommended for the researchers who are using fractals and 
multi-fractals in earth sciences including seismology. The papers included are 
important for both fundamental and applied research perspective. 


USA Donald L. Turcotte 
May 2015 


Preface 


1. Introduction 


The role of fractals in understanding the complex system in earth science has 
opened a new branch of science applied in geophysics, geology, geomorphology, 
environment, seismology, etc. (Mandelbrot 1982; Feder 1988; Barton and La Pointe 
1995; Turcotte 1997; Dimri 2000, 2005a, b; Dimri et al. 2012). Time series gen- 
erated by earth sciences observation are in general having long-term persistence 
with certain degree of correlation. Such correlation can be mapped in terms of 
fractal dimension for homogeneous models as the mono-fractals. For heterogeneous 
functions, mono-fractals are not applicable and multi-fractals have been proposed. 
A typical geological structure such as fault shows a relation with fractal in terms of 
fault length and displacement of fault zone thickness and throw. Methods such as 
Detrended Fluctuation Analysis (DFA) and Multi-fractal Detrended Fluctuation 
Analysis (MDFA) are central to many chapters applied to geophysical well-log 
data, seismological data and fire time series. 

In Chapter “Scaling Laws in Geophysics: Application to Potential Fields of 
Methods Based on the Laws of Self-Similarity and Homeogeneity”, Prof. Fedi sys- 
tematically described the interpretation of potential field data in a frequency domain. 
Initially, a relation between random sources (density/susceptibility) and their fields 
was established which led a spectral analysis method. This method was extensively 
used by geophysicists for the estimation of depth or thickness of sedimentary basin 
from gravity and magnetic data. Later, it was found that the source is not a random; 
rather it is a fractal distribution, and hence, a new method called scaling spectral 
analysis method was proposed. This method finds its application in interpreting 
gravity and magnetic data, e.g. Chapter “Scaling Laws in Geophysics: Application to 
Potential Fields of Methods Based on the Laws of Self-Similarity and Homeogeneity 
” of this book. On the other hand, the author has elegantly described the homogeneous 
source distributions called “one point” distributions. Dr. Fedi has proposed, as a 
multi-homogeneous model, having a variable homogeneity degree versus the posi- 
tion. He concluded that while mono-fractals or homogeneous functions are scaling 
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functions, multi-fractal and multi-homogeneous models are necessarily described 
within a multi-scale data set and recommended that specific techniques are needed to 
manage the information contained on the whole multi-scale data set. 

In another Chapter “Curie Depth Estimation from Aeromagnetic for Fractal 
Distribution of Sources”, Dr. Bansal et al. described scaling spectral method for 
estimation of Curie depth. The earth’s magnetic field is used to find depth of 
anomalous sources as shallow as few metres to tens of km. The deepest depths 
found from magnetic field sometimes correspond to Curie depth a depth in crust 
where magnetic minerals lose their magnetic field due to increase in temperature. 
Estimation of depth from magnetic/aeromagnetic data generally assumes random 
and uncorrelated distribution of magnetic sources equivalent to white noise distri- 
bution. The white noise distribution is assumed because of mathematical simplicity 
and non-availability of information about source distribution, whereas from many 
borehole studies around the world, it is found that magnetic sources and also other 
physical properties such as density, conductivity, etc. follow fractal distribution. 
The fractal distribution of sources found many applications in depth estimation 
from magnetic/aeromagnetic data. In this chapter, Curie depth estimation from 
aeromagnetic data for fractal distribution of sources has been presented. 

In Chapter “Fractal Faults: Implications in Seismic Interpretation and 
Geomodeling”, Dr. Ravi Prakash Srivastava has combined the concept arrived 
from two different words, one called “Fractals” that follows scaling law seen in 
many natural phenomenon and the other called “Faults” encountered in geological 
studies. Author presented most updated use of fractal concept in geological 
understanding of faults and their significance in geological modelling of hydro- 
carbon reservoirs. It is reported in many studies that the faults show fractal/scaling 
behavior in terms of fault length and displacement fault zone thickness versus fault 
throw. However, large faults are identified and integrated in seismic data, but same 
is not true when they are small and hence below the resolution limit of seismic data. 
Thus, author is able to map smaller faults with geological knowledge of the faults in 
the reservoir model studies using a concept of fractals. 

In Chapter “Detrended Fluctuation Analysis of Geophysical Well-Log Data”, 
Subhakar and Chandrasekhar used geophysical well-log data as it provides a unique 
description of the subsurface lithology as well as they represent the depositional 
history of the subsurface formations, vis-a-vis the variation of their physical 
properties as a function of depth. For this, DFA technique has been described and 
applied to gamma-ray log, sonic log and neutron porosity log of three different 
wells from west coast of India. The whole chapter is organized in seven sections 
dealing with geology of the area, data used in the present study, basic theory of 
DFA, methodology used, discussion and interpretation of results derived and finally 
the conclusion of the studies. 

Drs. Banerjee and Nimisha Vedanti in Chapter “Fractal Characters of Porous 
Media and Flow Analysis” proposed a fractal model of continuum percolation 
which quantitatively reproduces the flow path geometry. Porosity is an important 
property of geological formation and a complex function of many variables which 
control fluid flow. The variables include essentially the characteristics of pore 
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structures such as type, size, shape and arrangement of pores; pore space connec- 
tions; area of pores open for flow; and tortuosity of the flow paths and composition 
of pore. There is no general framework which explains the fluid flow through the 
fractured and complex subsurface geometry. In fact, a direct measurement of flow 
through complex permeable media is time taking; hence, an analytical model has 
been recommended. Therefore, the aim of this chapter is to develop a simple 
analytical model based on the medium structural characteristics to explain the flow 
in natural fractures. The authors showed that the pattern of fracture heterogeneity in 
reservoir scale of natural geological formations looks as the distributed self-similar 
tree structures. 

In Chapter “Estimation and Application of Fractal Differential Adjacency 
Segregation (F-DAS) Scores in Analysis of Scanning Electron Micro Graph (SEM) 
Imageries Towards Understanding the Adsorption unto Porous Solids”, Prof. Das 
et al. presented Fractal Differential Adjacency Segregation (F-DAS) which is dif- 
ferent from the conventional approach for estimation of fractal dimension 
(FD) using box counting method. The box counting method is one of first methods 
to estimate fractal dimension and is also simple to use. However, for heterogeneous 
system, mono-fractal seems to be inappropriate. In this chapter, the limitation of use 
of mono-FD is explained and thus extended the concept of Fractal dimensioning 
into lower scale segregation levels and evaluating their differential scores. In this 
approach, F-DAS scores are estimated for each of the image pixels of scanning 
electron microscopy (SEM) imageries using the arithmetic means of the grey levels 
of the adjacency pixels enclosed by the box used for counting in the conventional 
methods. The authors claimed that the present analysis provides better under- 
standing of variability of the system (in this case, adsorbents), unexplored by 
qualitative analysis of SEM imageries, as well as the functional groups using 
Fourier transform infrared spectroscopy. 

Dr. Padhy in Chapter “The Multi-fractal Scaling Behavior of Seismograms 
Based on the Detrended Fluctuation Analysis” has explained the use of the 
multi-fractal scaling properties of seismograms in order to quantify the complexity 
associated with high-frequency seismic signals and hence to characterize medium 
heterogeneities at different scales. He recommended the MDFA method is capable 
of characterizing the multi-fractality of earthquake records associated with fre- 
quency- and scale-dependent correlations of small and large fluctuations within 
seismogram. The multi-fractal nature of earthquake records has been explained by 
computing the generalized Hurst and mass exponent and multi-fractal singularity 
spectrum. One of the findings is that the degree of multi-fractality decreases with 
increasing frequency, and is generally more for the time period windowing domi- 
nant seismic phases in the seismogram. 

Finally in Chapter “Fractal Methods in the Investigation of the Time Dynamics 
of Fires: An Overview”, Prof. Telesca reviewed the fractal methods applied to fire 
point processes and satellite time-continuous signals that are sensitive to fire 
occurrences. Fires represent one of the most critical issues in the context of natural 
hazards. Most of these fires could be natural as well as anthropic. Sometimes, 
summer drought can influence the ignition and spread of devastating fire. 
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The author gave very good description of methods to find the temporal distri- 
bution of sequence of fire. These methods include coefficient of variation (CV), 
Detrended Fluctuation Analysis (DFA) and Multi-fractal Detrended Fluctuation 
Analysis (MDFA). Other techniques such as Fano Factor (FF), Allan Factor(AL) 
and Count-based Periodogram (PG) for fire sequence as the count process has also 
been described through fire satellite time series. Author also applied these methods 
to know the information about the “health” of vegetation which is very important 
and can be used to find the status of vegetation after the big forest fire. 
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Scaling Laws in Geophysics: Application 
to Potential Fields of Methods Based 

on the Laws of Self-similarity 

and Homogeneity 


Maurizio Fedi 


Abstract We analyse two classes of methods widely diffused in the geophysical 
community, especially for studying potential fields and their related source distri- 
butions. The first is that of the homogeneous fractals random models and the second 
is that of the homogeneous source distributions called “one-point” distributions. As 
a matter of fact both are depending on scaling laws, which are used worldwide in 
many scientific and economic disciplines. However, we point out that their appli- 
cation to potential fields is limited by the simplicity itself of the inherent 
assumptions on such source distributions. Multifractals are the models, which have 
been used in a much more general way to account for complex random source 
distributions of density or susceptibility. As regards the other class, a similar 
generalization is proposed here, as a multi-homogeneous model, having a variable 
homogeneity degree versus the position. While monofractals or homogeneous 
functions are scaling functions, that is they do not have a specific scale of interest, 
multi-fractal and multi-homogeneous models are necessarily described within a 
multiscale dataset and specific techniques are needed to manage the information 
contained on the whole multiscale dataset. 


1 Introduction 


The Earth is a heterogeneous medium, meaning it is characterized by a complex 
distribution of its physical properties. In general, a proper representation of its 
complexity has to explicitly refer to the specific scale of observation of a given 
phenomenon. To this regard, a special class of phenomena share a scaling property: 
that of not exhibiting any characteristic scale. The behavior of scaling phenomena 
assesses that all the scales within some range are equally important and the signals 
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at different scales are each one related to another. Different types of scaling are 
known in the literature, including the long-range dependency, self-similarity, fractal 
statistical similarity (monofractals) and local self-similarity: multifractals and infi- 
nite divided sequences. Since there is no characteristic scale, in many areas of 
science it has been useful integrating models across a broad range of scales. 

Earth’s physical parameters, such as magnetic susceptibility, resistivity or others, 
have been studied within the frame of fractals. At the same time, potential field 
anomalies (gravity, magnetic, self-potential) have been studied under the frame- 
work of scaling laws based on the homogeneity law. In particular, multiscale 
methods, such as the continuous wavelet transform (CWT), will be specifically 
designed to study field anomalies exhibiting such scaling properties. 

In this paper, we will try to analyse these two kinds of scaling laws in parallel, in 
order to define their similarities, limitations and advantages for determining the 
properties of the potential field sources. 


2 The Fractal Scaling-Law in Geophysics 


Random fractals are complex signals that, at first sight, seem not to deserve any 
useful information, due to their too strong erratic behavior. However, valuable 
information may be achieved if, instead of exploring the signal at a fixed scale, we 
try to search for a relationship existing between the observations of the signal at 
different scales. Scale-invariant fractals, in fact, are uniquely characterized by a 
single parameter, the fractal dimension, which gives us a description of its type and 
complexity. 

Consider a signal Y(x) and write the self-similarity equation (e.g. Mandelbrot 1983): 


Y(ax) =a"Y(x), a>0,x>0, (1) 


where Y(x) is a self-similar homogeneous fractal (or monofractal) with 0 < H < 1. 
The value of H, referred to as the Hurst parameter, will catch the degree of 
self-similarity of Y(x). For self-affine processes, the local properties are reflected in 
the global ones, resulting in the relationship: 


D+H=q+1 (2) 


where D is the fractal dimension and q is the space dimension (e.g. Russ 1994; 
Muniandy et al. 2003; Gneiting and Schlather 2004). 

Y(x) is indeed equal to its scaled version Y(ax) after the normalization by a", so 
that the phenomenon, observed at whatever scale, will present the same information 
and may be explained by the same model, no matter the scale. 

However, for more complex phenomena, reasonable models explaining a given 
phenomenon at a given scale may be different from those found at a different scale. 
In this case we may write (for instance, Arneodo 2000): 
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Y(ax) = a4 y(x), a>0,x>0 (3) 


where H(x) is now a function of x, so as to allow for more complex scaling 
behaviors. Consequently, the process is said to be multifractal and the function 
H(x) is referred to as Holder exponent. While self-similarity describes space scale 
invariance, which is the permanence of a structure on all scales ruled by a single 
law independent of x, multifractality captures different behaviors on small and large 
scales: characteristic features of the signal (for instance, bursts) are detectable on all 
scales, but they obey different laws depending on x. 

Real-world phenomena may be rarely described in terms of simple deterministic 
fractal models, but similarity can hold on several scales in a statistical sense, leading 
to the notion of random fractals. As an example, random fractals are the model 
assumed for well logs: Marsan and Bewan (1999) evidenced a multifractal distri- 
bution for the P-wave sonic velocities recorded at the German Continental Deep 
Drilling Programme (KTB) main borehole and for the gamma log. Fedi (2003) 
analysed susceptibility data from the deep KTB log and found that a multifractal 
model was more appropriate than a monofractal one for the statistical modelling of 
these signals. Power spectra or variogram analysis are in fact useful for the char- 
acterization of the signal up to the second-order statistics, but may fail in explaining 
more complex structures. Fedi et al. (2005) extended this kind of analysis to dif- 
ferent geophysical KTB well logs (density, magnetic susceptibility, self-potential 
and electrical resistivity) and found consistent information about the KTB well 
geological formations. They characterized the lithological changes of the drilled 
rocks and identified zones of macro- and microfractures by using a regularity 
analysis, which maps the measured logs to profiles of Holder exponents, or regu- 
larity. Regularity generalizes the degree of differentiability of a function from 
integer to real numbers and is useful to describe algebraic singularities related not 
only to the classical model of jump discontinuity, but to any other kind of ‘edge’ 
variations. 


3 Estimation of the Source Properties Using the Fractal 
or the Spector and Grant’s Scaling-Law 


An important issue is linking in a coherent way the complex behavior of the 
physical properties measured with well logs to other physical quantities, which are 
related to those parameters, but are independently measured. More clearly, we may 
define how the complexity of the medium, revealed from well logs of susceptibility, 
density or wave speed, is mapped to the complexity of fields, like the magnetic 
field, the gravity field or the seismic wavefield. 

Herrmann (1997) evidenced an inhomogeneous scaling for both well-log 
acoustic waves and reflectivity, so implying that the singularity structure is carried 
out from space to space-time. Other authors (for instance, Gregotski et al. 1991; 
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Pilkington and Todoeschuck 1993; Maus and Dimri 1994; Bansal et al. 2010) 
interpreted susceptibility and density logs in terms of scaling sources and consid- 
ered their magnetic or gravity fields as scaling quantities, with fractal dimensions 
related in a simple fashion to the fractal dimension of the source parameters. Naidu 
(1968), assuming a homogeneous random density sources confined to a thick 
semi-infinite medium, derived the relationships between the spectrum of the sources 
and the spectrum of the related gravity or magnetic fields. Pilkington and 
Todoeschuck (1993) assumed a 3D fractal magnetization, with power spectrum 


Sj(kxs ky, ke) = (ke + ky + k2)%/?, («<0O), and obtained for the field spectrum P: 


“Hay Bah ) 


P(s, hy) = 


where k,, ky, k are angular frequencies; s = (k2 fe ey" /2 is the radial wave number; 
h, is the depth of the half-space; P is the power spectrum of the magnetic field 
vertical component observed on the plane z = 0. 

Note that, by this approach, the scaling exponent is assumed to vary in the fractal 
range and, in turn, to depend on geology (Bansal et al. 2011). 

However, Fedi et al. (1997) have shown that other kind of correlated sources, 
such as the magnetically homogeneous blocks, inherently assumed by Spector and 
Grant (1970), also originate fields whose scaling exponent is in the fractal range, 
being it equal to about 3. To be clear, Spector and Grant (1970) modelled the 
susceptibility distribution using independent ensembles of rectangular, 
vertical-sided parallelepipeds. Each ensemble is characterized by a joint frequency 
distribution for depth (hf), width (a), length (b), thickness (7), magnetization 
intensity (j/ab) and direction cosines of magnetization. Under the hypothesis that 
the parameters vary independently and that they are randomly and uniformly dis- 
tributed about their mean values, the authors derived the well-known analytical 
expression of the radial power spectrum of the reduced to the pole magnetic field: 


P(s) = (2) “"PC(s,a,5)G(s,7)P(s,h) (5) 
where: 
7 =o 7 a = (a); b = (b); hain (6) 


sin(as cos 0) sin(bs sin 0) 


S(s, 0,a,b) = 
(s,0,4,b) as cos 0 bs sin 0 
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G(s,7) = (( - ey’) aie (9) 


As, h) = ((e™)’) = goa Sims) (10) 


6 = tan ‘u/v represents the angular direction of the radial frequency vector and the 
symbol <> stands for ensemble average. 

Now, the factor C was shown in Fedi et al. (1997, Fig. 1) to clearly follow a 
power law of exponent nearly equal to —3 for intermediate-to-large values of a, 
while for small values of a the spectrum is generally flat at low-to-intermediate 
wavenumbers and continues to follow the above power law only at large 
wavenumbers. So, especially for large thicknesses and when measurements are 
close to the source, the factors Q and G will both tend to 1 and the spectrum decay 
of the logarithm of power spectrum will be only determined by the power law of 
C (of exponent equal to —3), according to: 


P(s)=(2) "2C(s,a,B) (11) 


It is important to stress that it is not correct to consider the Spector and Grant 
model a random uniform uncorrelated distribution of sources (e.g. Pilkington and 
Todoeschuck 1993; Ravat et al. 2007; Bansal et al. 2011). As a matter of fact, this 
model is a statistical distribution of homogeneous sources (prisms) under the 
hypothesis that the source parameters (thickness, width, depth and others) vary 
independently and that they are randomly and uniformly distributed about their 
mean values. Hence, being the sources of the ensemble homogeneous, this is not 
yet an uncorrelated source distribution, but an ensemble of sources, which are, at 
least locally, correlated. The misunderstanding may be explained in this way: even 
though the magnetization factor j” is a constant (i.e. a white power spectrum), the 
magnetization distribution must instead be defined by the whole product 
F-C(s, a, b)T(s,t), which is a red power spectrum, as discussed above (Eq. 11). 

One more misunderstanding (e.g. Ravat et al. 2007; Bansal et al. 2011) is on 
using the definition of C(s,a,b) for making a correction factor (a ke spectral cor- 
rection) no matter the case. This is, in practice, a deconvolution of the spectrum for 
C, which was originally proposed in Fedi et al. (2007) and Quarta et al. (2000). 
However, these last authors made it well clear (Fig. 1 in Fedi et al. 2007) that the 
correction allows an improved estimation of h only when anomalies are big enough 
to justify, roughly, a homogeneous-block kind of source. Therefore, it should not be 
applied when @, b are too small or, in other words, when most anomalies appear to 
have in the map a too short wavelength. 

To explain a possible spatial variability of the scaling exponent, a relation to 
different geology units was invoked (Bansal et al. 2011). However, there is still a 
large amount of ambiguity, as shown by Quarta et al. (2000) who have studied the 
behavior of fields generated by sources other than fractal and have clearly shown 
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(b) 


Magnetic Field (nT) 


Magnetic Field (nT) 


Fig. 1 A fractal 3D homogeneous distribution of magnetic susceptibility (SI), with fractal 
dimension D = 3.5 (a) and its generated magnetic fields measured at: a negligible distance (b), 
0.5 km above the surface (c) and 2.5 km above the surface (d). As the altitude increases, the 
information about the fractal distribution is progressively lost 


that their spectral scaling exponents can assume values in the so-called fractal 
range. More specifically, sources with intermediate characteristics between Naidu’s 
(uncorrelated source distribution) and Spector and Grant’s (locally correlated source 
distribution) models will produce fields whose power spectra scaling exponent 
depends on the ratio between the horizontal extent of the source and the field 
sampling-interval. 

So, to verify the goodness of a fractal source model for experimental data, a 
simple test could be made, examining the properties of the whitened field at source 
level, since it must be random when any correlation is removed in the amplitude 
spectrum. 

More recently, Bouligand et al. (2009) have proposed a model of depth to the 
Curie temperature isotherm from magnetic anomalies in the western United States, 
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based on a fractal model. Their spectral model has a fractal exponent of 3.0, by 
which the depths to the bottom of magnetic sources were close to the depths to the 
Curie temperature isotherm derived from heat flow data. However, from the 
inherent spectral ambiguity discussed above, it is clear that there is no clear evi- 
dence that the source distribution is fractal, being it explained well also by a Spector 
and Grant’s model based on relatively homogeneous magnetized blocks. 


4 The Homogeneity Scaling-Law in Geophysics 


We have described above how scaling laws, of the fractal or multifractal type, are 
important to model the Earth as a scaling medium and to describe its degree of 
complexity. In particular, we have seen that either monofractal or multifractal 
models have been used to describe the well logs of susceptibility, seismic wave 
speed and other physical properties of the Earth’s rocks. 

In addition, we have mentioned in the previous section some studies linking, in 
the Fourier domain, geophysical quantities to physical rock parameters for fractal 
models. This shows that, besides modelling the direct information coming from the 
logs of physical properties of the rocks, geophysics may provide much wider and 
depth-related information through the study of measured quantities, which are 
related to such properties. 

In particular, spectral methods could allow the fractal properties of the source 
distribution to be estimated from the field. As already said, however, this is possible 
only when the field is studied very close to the source distribution, so that the 
Spector and Grant’s model becomes (Eq. 11): E(s)=(4) "PC (s, a,b). 

This is however also true for other models, such as that of the half-space fractal 
random noise in Eq. (4) (Naidu 1968; Pilkington and Todoeschuck 1993), which, in 
fact also assumes the form of a power law: 


8(—a — 1)! 


: m(—a)! 


II 


riot (12) 


In order to illustrate better this concern, let us use half-space fractal random noise in 
order to generate a fractal 3D homogeneous distribution, with fractal dimension 
D = 3.5 (Fig. la). The corresponding magnetic field (Fig. 1b), when measured at a 
negligible distance, is, in effect, very similar to the source distribution, showing an 
irregular behavior typical of the fractal model. It is so consistent with its power 
spectrum, fully reflecting the fractalness of the source distribution. 

However, when we are only a 0.5 km above the surface, the complexity of the 
field disappears(Fig. lc) and when we are at 2.5 km it is completely lost and the 
field has a nearly regular shape (Fig. 1d): we can only argue that the source volume 
contains magnetized sources, but any other information about the fractal distribu- 
tion can no longer be obtained. 
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This behavior is because the convolution factor related to the depth (Eqs. 4 and 
5) is clearly dominant and obscures the other one. So, being this last the most 
typical situation, we must resort to some other laws, better describing the source 
distribution. To this end, we will illustrate now the use of other types of scaling 
laws, widely used in applied geophysics, owning to the homogeneity law 


Ff (ax, ay, az) = a’f (x,y,z) a> 0, (13) 


where n is the degree of homogeneity and fis a homogeneous function (e.g. Florio 
and Fedi 2013). As for the fractal scaling law, a key role is again played by the 
mathematical properties related to the convolution. In fact, geophysical potential 
fields are mathematically expressed as the convolution of a source-density property 
with some specific operator. For instance, consider the total-magnetic-field: 


1 f\s@)| 
4nofar} [Ri 
Vv 


T (tg) = dV, (14) 


where J is the magnetization intensity, V is the source volume, R is the distance | 
r, — r| between the observation point, at ry, and the source points, at r in V, and 
f and t are, respectively, the main field and magnetization unit-vectors. Since the 
integral has itself the form of a convolution, its Fourier transform may be simply 
expressed in the harmonic region as (Maus and Dimri 1994; Caratori Tontini 
et al. 2009): 


~ _(f-k)(t-k); 
T= ae J (15) 


where k is the 3-D wave vector, and T,J are, respectively, the Fourier transforms of 
the magnetic field and of the magnetization. 

So, the convolution actually maps the source-density distribution complexity 
into that of the magnetic field, but this mapping is strongly dependent on the 
altitude, that is on the scale of observation. Obviously, similar considerations also 
occur for other integral equations, connecting, for instance, the mass-density dis- 
tribution to the gravity field or the electric conductibility to the electric and mag- 
netic fields. 

If one tries to classify physical source distributions in terms of their complexity, 
we have so far used statistical models of growing complexity, from fractal to 
multifractal models. On the other hand, we may use other classes of source dis- 
tributions to deal with a simplified medium and so describe the source distribution 
at different levels of complexity. These models may be generated solving the 
integral (14) or Eq. (15) for homogeneous-density distributions with very simple 
geometrical shapes, such as spheres, infinite cylindrical or planar distributions. 

Many of these simple source distributions fall in the class of so-called one-point 
sources (Fig. 2), which gained in the last two decades a strong popularity since their 
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Fig. 2 One-point source models. The magnetic fields of one-point source models obey the 
homogeneity relation (13) and are each one characterized by a homogeneity degree, n, which may 
be one of the integer values {—3, —2, —1, 0} 


fields obey the homogeneity relation (13). The one-point sources may be defined as 
source distribution being entirely characterized by the coordinates of only one special 
point. For instance, it is well known that the gravity field from a homogeneously 
dense sphere is equivalent to that of a mass distribution concentrated at its centre, 
provided the total mass is the same. Hence, the field of such sphere is completely 
defined by the coordinates of its centre, so that the main interpretative task will be 
finding a reasonable estimation of them. Other examples of one-point sources are the 
dipole for magnetic fields, infinitely extended pole and dipole lines, semi-infinitely 
extended tabular sources for modelling sheet or dykes, and the semi-infinite block 
model used to simulate the anomalies for contact-like geological structures. All these 
sources may be approximated by mass (or magnetization) distributions concentrated 
at their top or at their centre, depending on the source-shape. 

Note also that for such homogeneous sources, the estimated n will be one of the 
integer values {—3, —2, —1, 0} and we will be assured of a clear meaning for the 
estimated depth to the source: it can be referred to the source centre (sphere, 
horizontal cylinder, sill) or to the top (any other source). This feature makes 
somewhat different the fractal scaling-law (Eq. 4), in which the Hurst exponent may 
be any real value 0 < H < l(e.g. Turcotte 2011) and the homogeneity law (Eq. 12), 
in which the homogeneity degree is only integer and within the restricted set: {—3, 
—2, —1, 0}. However, it is important to stress that it has been recently shown that 
homogeneous sources and corresponding fields may exist, whose homogeneity 
degree may be fractional as well (Fedi et al. 2012; Fedi et al. 2015). They are 
characterized by having some intermediate properties between those corresponding 
to two subsequent one-point sources of integer n, as shown in Fig. 3, where the 
gravity anomalies are drawn for a pole (Fig. 3a), its fractional integration of order 
a =—0.5 (Fig. 3b) and for a bottomless vertical line mass (Fig. 3c). Basing on that, 
we may so extend the range of allowed homogeneity degrees to any real value 
within the range: —3 <n <0. 
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Fig. 3. The gravity anomalies due to a pole source (a), to a fractional order a = —0.5 integration of 
a pole source (b) and to a bottomless vertical line mass (ce) 


5 Estimation of the Source Properties Using 
the Homogeneity Scaling-Law 


In a previous section, we saw how estimating the source properties (Hurst or Holder 
exponent, depth to the source) using the scaling law related to the fractal or to the 
Spector and Grant’s model. Similarly, we now here briefly describe methods used 
in applied geophysics to estimate the homogeneity degree n and the depth to the 
sources thanks to the homogeneity scaling-law. It may be shown (e.g., Thompson 
1982) that the homogeneity relation (13) implies that homogeneous functions sat- 
isfy also the Euler differential homogeneity equation: 


of of of 


x %) + 50 yo) + oe zo) = —nf, (16) 


provided f is a differentiable function. This equation is particularly useful, since, 
differently than the homogeneity Eq. (13), it expresses the homogeneity properties 
of the fields with a differential equation or, in other words, in a local way at each 
point B(x, y, z) in the harmonic region. As such, it has been widely utilized using a 
local approach, i.e. estimating the unknown source parameters within a moving 
window centred at any point. This approach is referred to as Euler Deconvolution 
(Reid et al. 1990). Note that within the window the underlying hypothesis is that the 
homogeneity degree is assumed not to vary within the window, i.e. it is assumed 
that any source body giving rise to the sampled magnetic or gravity field in the field 
is simple enough to be represented by one singular point e.g., point mass or dipole, 
line source, sheet or dyke, and thick contact. 

To face the problem of the Euler differential Eq. (16) or of the homogeneity 
Eq. (13), different strategies have been adopted. We briefly group there some of the 
most used, in two categories, that of the monoscale methods and of the multiscale 
methods. 
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5.1 Monoscale Methods 


For monoscale methods, we mean the technique based on solving the Euler dif- 
ferential homogeneity equation with data only at the single level of the measure- 
ments. The first of them, shared by many Euler deconvolution algorithms (e.g. 
Thompson 1982; Reid et al. 1990; Barbosa et al. 1999; Mushayandebvu et al. 2001; 
Stavrev and Reid 2007; Fedi et al. 2009), mainly consists in the use of Eq. (4), 
slightly modified to account for an unknown background effect b: 


of 


Of 
Ox +3 


of 
Oy = 


ag zo) = —n(f +b). (17) 


(x — xo) yo) 


As we said, Euler deconvolution algorithms are based on solving a linear system 
of Euler’s Eqs. (17) for the unknowns Xo, yo, Zo, n and b, for a set of observation 
points covered by a moving window of size W, by assuming that the following 
approximated equation holds: 


f(r)~Fy(r) for any pointin W, (18) 


where F;,(r) is the best homogeneous field of integer degree n approximating the 
field f at points in the window W. The method allows for a different kind of 
homogeneous field Fy, at each window position, provided it is a good local 
approximation for f. 


5.2. Multiscale Methods 


Other methods deal instead with a multiscale analysis of potential fields, which are 
based on solving Eqs. (13, 16) with data at more altitudes. Within this category we 
find methods based on the CWT (e.g., Moreau et al. 1997; Sailhac and Gibert 2003; 
Saracco et al. 2004; Fedi et al. 2010). For potential field analysis, poissonian 
wavelets are typically used, making the method stable because the coefficients of 
the CWT are completely determined by the upward continuation of the field to the 
considered scales. Fedi and Cascone (2011) defined a more general approach, called 
Composite CWT, which let potential fields to be studied with wavelets other than 
poissonian, such as Morlet wavelet, Gaussian wavelets and any other. In such case, 
the coefficients at the different scales will be no longer related to the upward 
continuation, but the transformation remains very stable with any wavelet order. 
Other multiscale methods, such as the Depth from EXtreme Points (DEXP), 
SCALFUN (Florio et al. 2009) or the multiridge analysis (Fedi 2007; Fedi 
et al. 2009; Quarta 2009; Florio and Fedi 2014), make explicit reference to the 
upward continuation of any-order derivative of the field to the scales: / = z > Z,, 
where z,, is the altitude of the measurement surface. 
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For most multiscale methods the scale is equivalent to the altitude of continu- 
ation. A convenient tool to study the behavior of potential fields across the scales is 
the scaling function, defined by Fedi (2007) as: 


O logT Zz 
= =n i 
O logz Z— % 


t(z) (19) 


The scaling function is a dimensionless function of the altitude, which charac- 
terizes the scaling behavior of a homogeneous field (Fedi 2007). For homogeneous 
fields, the homogeneity degree n and the depth to the source z) may be estimated 
from t in a number of ways (Fedi and Florio 2006). For instance, altitudes and 
depths have an arbitrary zero level, so that they may be rescaled by a guess depth 
value ¢: 


Z=6 


Z- ZH 


(2) =n (20) 


and zp) may be estimated as the value of ¢ = zo making 7 the closest to a straight line 
with a zero slope. The homogeneity degree n may be instead estimated by the 
intercept of this line (see Fedi 2007, for further details). 

Euler Deconvolution is a local method, i.e. it characterizes the field punctually. 
As such, it may not only be applied to points on the measurement plane, but also to 
investigate the homogeneous scaling across scales, as recently proposed by Florio 
et al. (2014). Among multiscale methods the DEXP theory (Fedi 2007; Fedi and 
Abbas 2013) involves a nice property of homogeneous functions: considering a 
multiscale field f(z), i.e. the field at more scales (i.e., altitudes z; < z < z), and 
multiplying the data at each scale by z ””, the DEXP transformed field 


Q=f? (21) 


has as its extreme points (i.e. maxima and minima) occurring at the scale corre- 
sponding to the opposite of the depth to source Zo, that is at z = —Zp. 

Besides this, the DEXP transformation also has the property of yielding an 
image of the source distribution (Fedi and Pilkington 2012), which may be con- 
sidered as a preliminary step before more accurate inversion. Even though Q is not 
given directly in source-density units, it is very easy to it multiply by a physical 
constant allowing such units. Recent automatic algorithms to provide an automatic 
DEXP transformation are given in (Abbas and Fedi 2014; Abbas et al. 2014). One 
of the most striking features of these methods is that they provide an independent 
estimate of the homogeneity degree and of the depth to the source. The imaging 
properties of the DEXP transformation are shown in Fig. 4, where the Bouguer 
gravity anomalies of the Campanian area (a) are transformed in the DEXP 3D 
image (b), which well describes the source distribution in the underground source 
volume. The maxima of the DEXP image may be used to find the depth to the 
source position. 
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Fig. 4 The DEXP imaging. a Bouguer gravity anomalies of the Campanian area; b transformed 
DEXP 3D image. The source distribution in the underground source volume is well represented by 


the DEXP transformed image. The maxima of the DEXP image may be used to find the depth to 
the source position 


One of the main advantages of multiscale methods is their stability: they use in 
fact a composite upward continuation-derivative operator (Florio et al. 2009), which 
is designed to attenuate the noise, including the mutual interference effects from 
more sources, still maintaining a reasonable signal-to-noise ratio. The analysis is 
then performed along selected lines L in the 3D space, called ridges, which are built 
by joining the extreme points of the analysed field at different scales and assuming 
that the following approximated equation holds: 


f(r)~Fu(r) for points in a moving window W sliding along L (22) 


where again F;, is the best homogeneous field of integer degree n approximating the 
field in W. 
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6 Miultifractal Random Processes Vs. Inhomogeneous 
Fields 


Multifractal models allow managing random processes more complex than those 
related to a monofractal model and they so represent a useful generalization of the 
fractal theory. We will see now that the homogeneity scaling law can also be 
generalized, in order to deal with more complex source-density distributions than 
the one-point models. 

In fact, real-world source-density distributions are generally far from being well 
modelled by a single one-point source and n is expected to be a non-constant 
function of the spatial coordinates, being dependent on both the depth and kind of 
source distribution. In Euler Deconvolution, the observed scattering of solutions 
from subsequent windows is mainly due, even if not only, to approximate the 
measured field data with those from one-point sources. 

Steenland (1968) was probably the first to clearly illustrate that, for the most 
realistic cases of inhomogeneous potential fields, n is varying with the distance to 
the source, as shown in Fig. 5 for sources such as 3D and 2D thick prisms, thin 
plate 2D, thin circular disc. Note also that for such sources n not only varies but also 
assumes fractional values. This is completely different from one-point sources 
(solid lines: pole and infinite horizontal line) for which n is constant and assumes 
always an integer value. But this is no longer a problem, as we already said that 
homogeneous sources of fractional index might be defined in the framework of 
fractional potential fields (see Fig. 3). 

So, similarly to the generalization from monofractals to multifractals (Eqs. 1 and 
3), we need now to generalize the concept of homogeneous sources of integer 
homogeneity degree to that of multi-homogeneous sources of fractional and varying 
degrees. 

In practice, we have to replace the homogeneity Eq. (13) by: 


Ff (ax, ay,az) = a") F (x,y,z), (23) 


and the differential Euler homogeneity equation (Eq. 16) by: 


0 0 0, 
= PO w)+ F— 2) =—nlr— mf, (24) 


(x — xo) + 


Since n is now a varying space function, fis no longer homogeneous. In affinity with 
multifractals we can call it multi-homogeneous or, simply, inhomogeneous. 
Similarly, the homogeneity degree may be called multi-homogeneity degree. Fedi 
et al. (2013) proposed a method to define the source parameters of such inhomo- 
geneous fields analysing the scaling function of a multilevel dataset. They showed 
that the source parameters of inhomogeneous source distributions, such as the finite 
cylinder, might be successfully estimated owning to this new concept of generalized 
homogeneity. They showed that we may have access to the complex characterization 
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Fig. 5 The homogeneity degree of inhomogeneous fields versus altitude. For inhomogeneous 
sources, the homogeneity degree changes n versus the altitude and it is fractional. The curves 
relative to different cases are shown, namely the 3D and 2D thick prisms, the thin plate 2D, the thin 
circular disc. The only constant curves are related to one-point sources. The depth to the top is 
normalized in such a way that the half width is equal to 1 (from Steenland 1968, redrawn) 


of the sources by studying the inherent scaling properties of their fields across the 
scales z, < z < Z, where z, is the primary level of observation (or the lowest 
acceptable level of downward continuation) and z, is the upper acceptable level of 
upward continuation. 


7 Conclusions 


The key role of scaling properties in potential field theory is immediately clear if 
one considers the field of a complex source, such as a fractal, at different scales. On 
the other hand, the same may be said when observing the behavior of the fields 
related to another class of sources: the homogeneous-density or the uniform and 
homogeneous magnetization models owning to the one-point sources. 

In both the cases, as it is well known in potential field theory, potential or its 
derivatives are convolutions between source density and specific operators and 


16 M. Fedi 


actually map the source-density distribution complexity into that of the field. For 
fractal source distributions, due to the convolution properties and to the form of the 
convolution kernel, at increasing altitudes, such field tends to lose the irregular 
behavior and tends to be much more close to that of a relatively homogeneous-density 
source distribution, so as to be better interpreted by the models based on the 
multi-homogeneity theory, i.e. by a behavior similar to the black lines in Fig. (1). 

If we consider the asymptotic expansion of the potential field (e.g. Tikhonov and 
Samarski 2013), multipolar terms will change their role as the scale changes, 
obviously implying that the field scaling characterization will change accordingly: 
close to the source the scale invariance typical of the fractal source distribution will 
be observed, while far from the sources the multipole terms other than pole are very 
low and the potential will show the scale invariance typical of one-point spherical 
sources: n = —2 for the gravity field of a homogeneous-density sphere and n = —3 
for the magnetic field of a homogeneous and uniformly magnetized sphere. 

We also described a fourth class of source distributions, that is of sources whose 
complexity is such that we cannot describe it either as that of random source models 
(monofractals and multifractals) or as that of one-point functions. These fields are 
inhomogeneous, meaning that the homogeneity degree varies with the scale, and are 
very common in real world. Some examples are the source distributions due to 
complex but not yet random source distributions, as the field of a polyhedral or that 
of any irregular geometry. For these sources, scale invariance could be approxi- 
mately assessed by methods such as Euler Deconvolution, but only within a very 
restricted range of scales. One may however find true evidence for scaling in the 
asymptotic regions, i.e. within only the very large or the very small scales. 

Our main conclusion is that new methods are requested for handling with 
inhomogeneous fields, being the classical homogeneity law (Eq. 13) not applicable, 
unless than in an approximate way. In fact, similar to what happens for random 
source models, where the fractal scaling-law (Eq. 1) has been generalized into a 
multifractal law (Eq. 3), the homogeneity law must be replaced by a 
multi-homogeneity law (Eq. 23). Along this line, one method was recently pro- 
posed by Fedi et al. 2013; Fedi et al. 2015, proving to be successful when the 
interpretation of a multiscale dataset is developed in the framework of homoge- 
neous field models of fractional-degree. 
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Curie Depth Estimation 
from Aeromagnetic for Fractal 
Distribution of Sources 


A.R. Bansal, V.P. Dimri, Raj Kumar and S.P. Anand 


Abstract The earth’s magnetic field is used to find the depth of anomalous sources 
as shallow as few metres to tens of kilometres. The deepest depths found from the 
magnetic field sometimes correspond to Curie depth, a depth in the crust where 
magnetic minerals lose their magnetic field due to increase in temperature. 
Estimation of depth from magnetic/aeromagnetic data generally assumes random 
and uncorrelated distribution of magnetic sources equivalent to white noise distri- 
bution. The white noise distribution is assumed because of mathematical simplicity 
and non-availability of information about source distribution, whereas from many 
borehole studies it is found that magnetic sources follow random and fractal dis- 
tribution. The fractal distribution of sources found many applications in depth 
estimation from magnetic/aeromagnetic data. In this chapter Curie depth estimation 
from aeromagnetic data for fractal distribution of sources will be presented. 


1 Introduction 


Curie depth is a depth in the earth’s crust where ferromagnetic mineral changes to 
paramagnetic due to increase in temperature and generally no detectable magnetic 
field is observed below this depth in the crust. This depth may be very deep or 
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shallow depending on heat flow in the region and composition of rocks. The Curie 
depth in a region serves as a proxy for heat flow. Direct measurement of heat flow is 
carried out from bore holes in land and drilling cores from the ocean. The direct 
measurement of heat flow is very sparse because of higher costs for drilling, 
whereas proxy estimation from aeromagnetic/magnetic data provides homogeneous 
coverage. The magnetic surveys are easy to carry out and cheap in terms of cost 
involved in field surveys. 

The estimation of depth from magnetic data is an ambiguous process and there 
are many methods for its estimation. None of the methods provide reliable depth 
estimations. These methods are applied in space and frequency domain. The fre- 
quency domain methods are generally preferred as compared to the space domain 
(Gundmundsson 1966, 1967; Heirtzler and Le Pichon 1965; Neidell 1966; Naidu 
1968, 1970; Bhattacharyya 1967; Spector and Grant 1970; Treitel et al. 1971; Negi 
et al. 1986; Dimri 1992) because convolution operator changes to multiplication. 
The Fourier domain methods became very popular in depth estimation because of 
simplicity. The classical methods of depth estimation assume random and uncor- 
related distribution of sources which is equivalent to white noise distribution. 

From many boreholes it is found that source distribution follows random and 
fractal distribution of sources which is known as scaling distribution (Pilkington and 
Todoeschuck 1990, 1993, 2004; Pilkington et al. 1994; Maus and Dimri 1995a, b; 
Leonardi and Kumpel 1998; Bansal et al. 2010; Bansal and Dimri 2010). The power 
spectrum of scaling distribution is frequency dependent in contrast of white noise 
which is frequency independent and mathematically it is defined as: 


P(k) = Ak? (1) 


where P is the—power spectrum, k is the wavenumber, f is the scaling factor and 
A is constant. The values of scaling exponents represent degree of correlation and 
larger the value stronger is long range correlation. 

The values of scaling exponents due to source and field are related and estimated 
easily if one is known and vice versa (Maus and Dimri 1994). A known value of 
scaling exponents can be easily represented in different dimensions using simple 
formula (Maus and Dimri 1994). The # values are found to vary with space 
(Pilkington and Todoeschuck 1993; Bansal et al. 2010; Bansal and Dimri 2014) but 
still the exact relation with tectonic and rock formations is not yet established 
mainly due to limited studies. The use of scaling distribution of sources provides 
better depth estimation (Pilkington and Todoeschuck 1993; Maus and Dimri 1994, 
1996; Fedi et al. 1997; Bansal and Dimri 1999, 2001; Bansal et al. 2006a, b; Dimri 
2000; Dimri et al. 2003). The values of scaling exponents can be converted to 
fractal dimension using a simple formula (Mandelbrot 1982; Turcottee 1997; 
Bansal and Dimri 2005a). Fedi et al. (1997) have shown inherent power law due to 
Spector and Grant ensemble. 
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Curie depth estimations are carried out worldwide from aeromagnetic data using 
Fourier domain methods mostly assuming random and uncorrelated distribution of 
sources (Okobu et al. 1985; Tanaka et al. 1999; Chiozzi et al. 2005; Trifonova et al. 
2009). Few recent studies have claimed better depth estimation for fractal distri- 
bution of sources (Maus et al. 1997; Bouligand et al. 2009; Bansal et al. 2011; 
Salem et al. 2014). 


2 Conventional Centroid Depth Method 


In centroid method (Bhattacharyya and Leu 1975; Okobu et al. 1985; Tanaka et al. 
1999) Curie depth is estimated in two steps: (1) top depth of anomalous body and 
(2) centroid depths are computed from the power spectrum of magnetic field data 
and then these depth values are converted to Curie depths. Spector and Grant (1970) 
proposed a method to estimate top depth of assemblage of magnetic sources. In this 
method, power spectrum of total magnetic field is represented in terms of top depth 
and thickness of magnetic body (Blakely 1995): 


2 
P( kes ky) = 40 C2 (kes by) Oml?| Oe]? x (1 eH@Ra))" (2) 


where k, and k, are the wavenumbers in the x- and y-directions; C,, is a constant of 
proportionality; g,, is the power spectrum of the magnetization; ©,, and ©, are the 
directional factors related to the magnetization and geomagnetic field, respectively; 
Z, and Z, are the top and bottom depths of the magnetic sources. 

It is common practice in geophysics for converting 2-D power spectrum to 1-D 
by taking radial average. In this case terms @,, and @; become constant and gy is 
constant for random and uncorrelated distribution of sources. In case of radial 
averaging of power spectrum, and random and uncorrelated distribution of sources, 
Eq. 2 can be written as: 


2 
P(k) = Aye2& (1 7 ev Hara) | (3) 


where A, is a constant and for very thick magnetic body the right-hand side of Eq. 3 
contains only top depth and Eq. 3 reduces as: 


P(k) = AyeW 2M (4) 


Equation (4) is frequently used for finding the top depth of anomalous magnetic 
bodies. 
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The centroid depth of magnetic body is given as (Bhattacharyya and Leu 1975, 
1977; Okobo et al. 1985; Tanaka et al. 1999): 


The Curie depth is finally computed from centroid and top depth as: 


Zn = 225 -—Z, (6) 


The centroid method has become very popular for estimating Curie depth from 
aeromagnetic data and is applied to aeromagnetic data of many parts of world 
(Bhattacharyya and Leu 1975; Okubo et al. 1985; Tanaka et al. 1999; Okubo and 
Matsunaga 1994; Chiozzi et al. 2005; Dolmaz et al. 2005; Trifonova et al. 2009 etc.). 


3 Fractal Based Methods of Curie Depth Estimation 


Fourier domain methods became very popular in Curie depth estimation from 
aeromagnetic data because of their simplicity. However, these methods provide 
overestimation of depth values due to the assumption of random and uncorrelated 
distribution of sources (Pilkington and Todoechuck 1993; Maus and Dimri 1994, 
1996; Fedi et al. 1997; Bansal and Dimri 1999, 2001, 2005b, 2014). These methods 
are modified for scaling distribution of sources. The scaling distribution of sources 
is equivalent to fractal distribution of sources and these modified methods are called 
fractal based methods of depth estimation. 

Maus et al. (1997) proposed a method for estimating Curie depth for fractal 
distribution of sources where scaling exponents, top depth and thickness of mag- 
netic body are estimated simultaneously from power spectrum of magnetic field. 
Radial average of power spectrum is expressed in terms of scaling exponent and 
depth component as (Maus et al. 1997): 


2\ —1-f/2 
P(k) = C — 2kz, — tk — Bin(k) + In J [cos h( (tk) — cos(tw mi(1+ 5) ov) 
0 


where k is the wavenumber, z, is the top depth, t is the thickness of slab and £ is the 
scaling exponent due to source distribution, w is the wavenumber in vertical plane. 
Maus et al. (1997) found a value of f equal to 4 from aeromagnetic data of South 
Africa and Central Asia and Curie depths vary from 15 to 20 km. Bouligand et al. 
(2009) derived an analytical solution for solving Eq. 7 and found difficulty in 
simultaneous estimation of top depth and scaling exponents. Therefore, Bouligand 
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et al. (2009) fixed the value of scaling exponents for estimating top depth based on 
the shape of power spectrum of aeromagnetic data. Manual checking of estimated 
parameter is essential for a reliable estimation. Bansal et al. (2011) proposed a 
modified centroid method for the estimation of Curie depth from aeromagnetic data 
for scaling distribution of sources. This method computes Curie depth in two steps 
similar to classical centroid method for scaling distribution of sources. The top and 
centroid depths are computed by correcting power spectrum for scaling distribution 
of sources as: 


Top depth: 
In(kPP(k)) = Az — 2kz, (8) 
and Centroid depth: 
P(k 
In (x =) = A; — 2kzo (9) 


Bansal et al. (2011) also pointed out the difficulty in estimating scaling exponent 
and depth values simultaneous from inversion method and they fixed the value of 
scaling exponent equal to | corresponds to 1/f noises found from seismic velocity 
fluctuations and fault structures (Holliger 1996). Salem et al. (2014) also corrected 
their power spectrum before computing the top depth from the power spectrum of 
aeromagnetic data and applied to the magnetic data of central Red Sea. Their values 
of scaling exponent vary between 0 and 1.7 with an average value of 0.85 very 
close to 1 used by Bansal et al. (2011) in the estimation of Curie depth. The Curie 
depths estimated by classical and modified centroid method for fractal dimension 
sometimes have a large difference (Table 1). This method is successfully applied to 
the German, Indian and Nigerian aeromagnetic data (Bansal et al. 2011, 2013; 
Nwankwo 2015). 


Table 1 Comparison of depth values computed using fractal (scaling) and non-fractal 
(conventional) distribution of sources 


Block no. Estimated Curie depths % Higher 
Scaling Conventional 
33 47 30 

3 28 42 33 

21 35 51 31 

31 37 63 41 


The difference in the computed values may be more than 40 % in some of cases (Table 1). The 
difference in computed values also depends on values of scaling exponents used. Higher the value 
of scaling exponent used larger is difference in computed values 
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4 Case Study: Application of Modified Centroid Method 
for Fractal Distribution of Sources to Central India 


Bansal et al. (2013) carried out a detailed study of the Curie depth estimation in 
central India. The central India is tectonically complex and covers many geological 
entities, e.g. Deccan volcanic province, central Indian Tectonic zone, Godavari and 
Mahanadi intracratonic failed rifts, Chhattisgarh basin and the Proterozoic Eastern 
Ghat Mobile Belt. 

Aeromagnetic data over central India is compiled by the Indian Institute of 
Geomagnetism over a common elevation of 1.5 km (Rajaram and Anand 2003; 
Rajaram et al. 2009). We selected four blocks of dimension 200 km =< 200 km 
covering Eastern Dharwar, Godavari-Graben and Baster Craton. Centres of selected 
blocks are lying in Godavari-Graben (Blocks 2, 21) and Baster Craton (Blocks 3 
and 31), whereas a large size of block covers surrounding geological entities. The 
Curie depths are computed using conventional (Figs. | and 2) and centroid method 
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Fig. 1 Estimation of Curie depth for central India using conventional centroid method for block 2 
and 21. The /eft and right panels indicate the estimation of centroid and top depth, respectively 
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Fig. 2. Estimation of Curie depth for central India using conventional centroid method for block 3 
and 31. The /eft and right panels indicate estimation of centroid and top depth, respectively 


for fractal distribution of sources (Figs. 3 and 4). Godavari-Graben is a passive rift 
orthogonal to the east cost of India and the Curie depth is found to vary between 33 
and 35 km (Fig. 3). Deep seismic study has shown Moho depth of 37 and 42 km in 
the Eastern Dharwar Craton (Reddy et al. 2005) and Godavari-Grabens (Kaila et al. 
1990). Godavari-Graben region is found to have underplating of high density rocks 
(Behera et al. 2004; Rao 2002) and evolved during permo-carboniferous rifting. 
The region has undergone a number of volcanism in the past and these thermal 
episodes have a large effect on the Curie depth in Godavari-Graben. The centres of 
Blocks 3 and 31 are on the south part of Baster Craton and magnetic data also 
covers Godavari-Graben and Eastern Ghat mobile belts. The Baster Craton is 
bounded by Godavari-Graben in the west and Eastern Ghat mobile belt in the east. 
The Baster Craton contains vast traces of granites and gneisses with basement of 
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Fig. 3. Estimation of Curie depth for central India using modified centroid method for fractal 
distribution of sources for block 2 and 21. The left and right panels indicate estimation of centroid 
and top depth, respectively 


mainly Archean to mid-proterozoic. The Moho depth in the North of Baster Craton 
is found as 48 km from deep seismic studies (Mandal et al. 2013). In the south of 
Baster Craton, values of Curie depth are found to vary from 28 to 37 km (Fig. 4). 
The region of the lowest Curie depth includes Godavari-Graben and Eastern Ghat 
mobile belt. The south-west part of Eastern Ghat mobile belt is found to have lower 
values of Curie depth 26-27 km from an earlier study (Bansal et al. 2013). The 
Curie depths estimated in central India using fractal distribution of sources are 
found to be lower than the values computed using conventional method (Fig. 5, 
Table 1) and reasonably well while considering other tectonic and geophysical 
constraints. 
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Fig. 4 Estimation of Curie depth for central India using modified centroid method for fractal 
distribution of sources for block 3 and 31. The left and right panels indicate estimation of centroid 
and top depth, respectively 


5 Discussion and Conclusion 


Depth estimation from magnetic data is ambiguous due to Green’s equivalent layer 
problem. The information about source distribution is limited mainly due to limited 
information below the surface of earth from deep boreholes. Assumption of random 
and uncorrelated distribution of sources resulted in a simple relation between power 
spectrum and depth of magnetic sources. This method is commonly used for finding 
Curie depth from aeromagnetic data. The magnetic susceptibility distribution is 
found to follow random and fractal distribution. The fractal distribution of sources 
is incorporated in estimating depth from magnetic data. Many authors have claimed 
that methods based on fractal distribution of sources provide better depth estimation 
as compared to depth estimation based on white noise distribution (Maus et al. 
1997; Bouligand et al. 2009; Bansal et al. 2011). Fractal based approach suffers the 
limitation of simultaneous estimation of depth and scaling exponents. At present, 
prior fixation of one of the parameters provides a better estimation of the other 
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Fig. 5 The comparison of depth values estimated using conventional and modified centroid 
method 


parameter. There is no doubt about better depth estimations using fractal distri- 
bution of sources and many researchers are working on simultaneous estimation of 
scaling exponent and depth values. 
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Fractal Faults: Implications in Seismic 
Interpretation and Geomodelling 


Ravi Prakash Srivastava 


Abstract Nature is not random as we often assume for the simplicity of mathe- 
matical calculations. Scaling or power laws, also known as fractal behavior, are 
ubiquitous in nature, and the analysis of many physical properties of the earth 
shows fractal behavior. In this chapter attempt is made to integrate the geological 
understanding and fractal behavior of the faults to use this knowledge in practice. 
This understanding could potentially help to reduce the uncertainty and risk in the 
fault modelling and some of the properties associated with the faults such as 
transmissibility and shale gauge ratio. The study aims at understanding faults in 
hydrocarbon reservoirs, however, the concepts are universally valid and could be 
useful for the readers interested in seismology or mining-related studies as well. 


1 Introduction 


In this chapter, I am tasked to combine the understanding from two different worlds, 
one called ‘Fractals’ that follows though simple but a fascinating mathematics of 
self-similarity and power law observed in many natural phenomena and the other 
called ‘Faults’ encountered in geological studies. This chapter is by no means a new 
research finding, rather an account of the most updated use of fractal concept in 
geological understanding of faults and their significance in geological modelling of 
hydrocarbon reservoirs. 

The study about the mapping and modelling faults in Earth science is one of the 
most important issues, particularly because of its serious impact on the work pro- 
cess. For example, it can have severe impacts on the mining industry and can lead 
to fatal consequences if faults are not understood before commencement of the 
mining work. Unexpected faulting can cause dilution and ore losses in underground 
metal mines, shutdowns and delays in production in underground coal mines with 


R.P. Srivastava (DX) 
National Geophysical Research Institute, Uppal Road, Hyderabad 500007, Telangana, India 
e-mail: ravi.ngri@ gmail.com 


© Springer International Publishing Switzerland 2016 33 
V.P. Dimri (ed.), Fractal Solutions for Understanding Complex Systems 

in Earth Sciences, Springer Earth System Sciences, 

DOI 10.1007/978-3-319-24675-8_3 


34 R.P. Srivastava 


consequent severe financial losses, geotechnical hazards that impact upon safety, 
etc. (Dimitrakopoulos and Li 2001). 

In oil and gas industry, again faults are very important. In this case both, large 
and very small scale faults influence the reservoir model. Here, the faults clearly 
visible in seismic data are considered large, and many faults, which are extracted 
based on some seismic attributes (spatial correlation techniques, viz. semblance, 
etc.) but not resolved by the reflection seismic data, are considered small scale 
faults. The presence of faults does influence the fluid flow properties in a reservoir. 
The magnitude of influence is mainly dependent on the intensity of faulting and the 
relative permeability contrast of the fault permeability compared with the perme- 
ability of the undeformed reservoir. Many factors control the influence of the faults 
on fluid flow, some of them are: the degree of faulting and strain localization, the 
geometry of the faults, the lithologies that have been involved in the faulting, their 
mechanical properties during faulting and the diagenetic alteration of the fault rock 
during burial. 

Generally, for siliciclastic sediments the fault rocks have lower permeability than 
the undeformed sediments. It is crucial to model the faults’ impact on the fluid flow 
within the reservoir to achieve realistic fluid flow and drainage patterns within 
dynamic reservoir models. 

The motivation to use fractal theory in understanding the fault geometry comes 
from the field observations and has also been reported in many studies that the 
faults show fractal/scaling behavior in terms of fault length and displacement 
(Spyropoulos et al. 2002; Scholz et al. 1993), fault zone thickness versus fault 
throw (Childs et al. 2009). Also, it is important to mention that it is often easy to 
identify and interpret large faults in seismic data, but the interpretation and iden- 
tification of faults is not very easy when they are small and hence below the 
resolution limit of seismic data. Thus, an attempt is made to integrate fractal 
behavior of faults with geological knowledge of the faults to better model the 
smaller faults in the reservoir model (called geo-model). 

Some researchers have used simulation methods for fault populations based on 
various approaches including fractals which have been developed in the modelling 
of petroleum reservoirs (e.g. Gauthier and Lake 1993; Munthe et al. 1993; Chilés 
et al. 2000; Mostad and Gjerde 2000; Holden et al. 2003). The modelling results are 
often crosschecked with the well data if available. 


2 What Is ‘Fractal’, ‘Fault’ 


Assuming that the readers of this book may be from different disciplines of science, 
the two main terms used in this chapter ‘fractal’ and ‘faults’ are defined below for 
the sake of completeness. 
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3 Fractal 


The main theme of this volume deals with the fractal; hence, instead of repeating 
the definitions of fractals, an attempt is made to describe what does fractal mean 
when we talk about the faults in geology. Traditionally fractal is defined as a 
geometrical shape having self-similarity independent of scale. Often, we see that 
faults do show self-similar pattern, though not on all the scales. In addition, it is 
known fact that none of the natural fractals exhibit self-similarity on all scales 
except synthetic fractals. In this context, we are more interested in statistical defi- 
nition of fractal behavior of faults. Statistically, if certain property (measured 
quantity) of an object follows power-law relation or long-term dependence on 
another property of the same object, then we say that it is fractal. For example, it is 
known that there is a long-term relation between magnitude and frequency of the 
earthquakes (Dimri 2000, 2005a, b; Scholz 1997; Turcottee 1992). Here, both 
magnitude and frequency are the properties related to earthquake. Similarly, an 
attempt is made to explore which properties of the faults are fractal. Once, we 
succeed in establishing a fractal behavior, we can use it in statistical analysis e.g. in 
forecasting the number of faults, or interpolation or extrapolation of one property of 
the fault based on another observed property. 


4 Fault 


A fault is the result of shear failure within a rock mass causing a displacement of 
subplanar zone due to sliding motion. There are many types of faults, but a common 
and easy-to-understand fault is normal fault. The geometry of normal fault is shown 
in Fig. 1, with the annotation of commonly used terms to describe a fault. In a 
normal fault, the hanging wall moves downwards, while the footwall moves 
upwards. In a reverse or thrust fault, the hanging wall moves upwards, while the 
footwall moves downwards (Fig. 2). In a strike slip fault, the sense of motion on the 
fault plane is horizontal. Frictional processes within the fault zone produce sheared 
material known as fault rock, which is characterized by the different petrophysical 
properties than the surrounding rocks. 

Note that in case of normal fault, it is extension of the rock mass, whereas in case 
of reverse fault it is shrinking of rock mass. To be mathematical, one can say heave 
(h) is positive in case of normal fault and it is negative in case of reverse fault. 
However, geologists do not define it in this way. 

In Fig. 3 a schematic fault is shown to describe some common terms used in 
fault geometry and Fig. 4 shows natural fault found in nature. The principal stresses 
involved in different kind of faults are shown in Fig. 5. 
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Fig. 1 Normal fault, h heave, t throw, d displacement 


Fig. 2 Reverse fault, h heave, t throw, d displacement 


fe 
- ault core 


Damage zone, 
brittle part 


Damage zone, 
ductile part 


Fig. 3. Figure explains the common terms used to describe a fault geometry, modified from 
Fossen and Gabrielsen (2005) 
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Fig. 4 Natural fault, notice the brown shale smeared along the fault plane. Thin red line highlights 
fault plane. A simplified interpretation of the natural fault is shown in (b) which usually geologists 
produce and that goes into geo-model. Important is to notice the smear of the shale material along 
the fault plane (shown by arrow), which makes this fault sealing, and fluid flow is not possible 
across the fault in such cases 


5 Fractures 


Fractures are often associated with the faults, however, they occur independent of 
fault presence. Fractures also show similar characteristic as faults, and they have 
high significance in reservoir modelling. 

Fractures can be defined as mechanical breaks in rocks involving discontinuities 
or displacement across surfaces or narrow zones’. A conceptual model of fracture 
formation is shown in Fig. 6. 

There are three major types of fractures (Fig. 7) such as: 


1. Dilatant fracture including joints 
2. Shear fractures/faults 
3. Closing fractures/stylolite 


As the state of stress in the earth’s crust is compressive, dilatant fractures usually 
require elevated pore pressure in order to form in the sub-surface. Closing fractures 
tend to form more easily in carbonates than in siliciclastic rocks and become more 
prevalent at depth. Whereas faults can occur at any depth in the brittle crust. Most 
faults initiate and propagate as dilatant fractures, but not all fractures will 
become faults. Those interested in more details about the fractures may refer to 
(Gillespie et al. 1993). 
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Fig. 5 Fault types and governing stresses involved in the faulting process (modified from Ramsay 
and Huber 1987). Sigma J is maximum principal stress direction, sigma 2 is intermediate principal 
stress and sigma 3 is minimum principal stress direction 


In reservoir modelling (mostly in case of oil and gas reservoirs, with high 
economics involved), the term ‘fractured reservoir’ is used to describe the reser- 
voirs in which fractures enhance permeability (Dimri et al. 2012). This permeability 
enhancement is most important in formations with low permeability rocks, such as 
carbonates and tight sands. Fractured reservoirs often contain all the three kinds of 
fractures as described above, but dilatant fractures and faults are most important for 
the fluid flow in the reservoir. In most fractured reservoirs, although the fractures 
provide essential additional permeability to the reservoir, their contribution to 
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Fig. 6 Fracture formation process (/) undeformed rock mass, (2) fracture development due to 
extension of upper surface, simultaneously compression takes place in lower surface as shown in 
(3), later depending on the rigidity of the rock mass, fractured parts may fail and create several 
small faults due to vertical shear as shown in step (4) 


fracture porosity is small. Typically, in a producing reservoir, the fluids are drawn 
from rock matrix into the high permeability pathways provided by the fracture 


system. 
2 


Fig. 7 Types of fractures (/) Dilating fracture (2) and (3) are shear fractures or faults 
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6 Are Faults Fractal? 


There is fairly good amount of published research work which dwells on the fractal 
geometry of faults particularly in light of seismological studies (Sunmonu and 
Dimri 2000; Sengupta et al. 2011 and references therein). It has been shown that the 
faults and lineaments show fractal behavior. A closer look at the earthquake 
foreshock and aftershock analysis reveals that earthquake’s magnitude and fre- 
quency follow fractal behavior (Dimri et al. 2005).Extending the notion of earth- 
quake analysis to the faults, it is likely that the faults should show a power-law size 
distribution and a linear scaling of the strain with the length of the fault. 

In some instances, besides power-law distribution of the faults, some fault 
populations exhibit exponential laws for the distribution of the fault size. 
Spyropoulos et al. (2002) suggest that the power-law distribution and exponential 
distribution of the faults are the transition regimes between the two end member 
stages. These end member stages are an initial uncracked state and the saturated 
state in which the cracks are evenly spaced. According to the experimental and 
numerical modelling done by Spyropoulus et al. (2002) the power-law distribution 
appears only in the process of the undisturbed growth. In the processes of nucle- 
ation and of coalescence and saturation they observed an exponential law. Readers 
are encouraged to look at the paper by Spyropoulus et al. (2002) for the definition 
of nucleation, coalescence and saturation of the cracks in their experiment. 

Davy et al. (1990) based on the laboratory experiment mimicking India—Asia 
collision concluded that the continental faults are fractal in nature and explained the 
occurrence of large stable basins based on this assumption. Further, Endres et al. 
(2008) analysed several fracture attributes such as orientation, density or spatial 
distribution and the fracture length in greater detail. They proposed that, if the ratio 
of ‘long’ and ‘short’ fractures is constant over different scales, then it follows the 
law of scale invariance, implying that the relationship is fractal. This means that 
there is a power-law relation between those ‘long’ and ‘short’ fractures. In general, 
the theory suggests that various fault parameters are invariant with respect to scale 
or are ‘self-similar’, providing a model that can be used for predictive purposes. 
Based on empirical data analysis, fault size distributions (throw or length) can be 
described by a power-law model over a wide range of fault sizes such that: 


log(N) = k — Blog(L) (1) 


where: 

Nis the cumulative number of faults of size greater than or equal to fault size “L’. 

L is either length or throw of the fault. 

k is a function of the fault density and when k is high, the fault density is also 
high. 

fis the fractal dimension of the fault population that defines the relative number 
of large and small faults; when f is high, the number of small faults is high 
relative to the number of large faults. 
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Fig. 8 Cumulative number of 

faults versus fault length, 

plotted on log-linear scale 100 
(modified from Gupta and 

Scholz 2000) 
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Cumulative number of faults 


Fault length (km) 


Apart from the data obtained from the experimental study, many authors 
attempted to find scaling laws for natural faults. It is interesting to note that the 
researchers, who used data from only one set of fractures, reached distributions very 
similar to the results of Spyropoulos’ (2002) experimental study. The closest 
agreement seems to occur with the results of Gupta and Scholz (2000), who 
examined data from the natural faults in the Afar depression. They observed a strain 
regime transition from a power law to an exponential law for the frequency-size 
distribution (Fig. 8), when faults reach a certain density. For the faults of the Afar 
depression, this density equals 0.6 km of fault length per square kilometre. This 
corresponds with a strain of 6-8 %. Above this density, faults grow mostly by 
coalescence until the area is saturated. 

Thus, we believe that for all practical purposes it is safe to go with the 
assumption that the properties of the faults are fractal. 


7 Data Sampling 


As defined earlier, fractals are a statistical behavior of the data that follows power 
law. The problem in analysing the scaling laws of faults is spare sampling of the 
data to establish a reliable distribution law. Often, to overcome inadequate data, 
many researchers used data from different geological and tectonic settings, the 
distribution and scaling laws established with such data may depend on different 
material properties. For a confident analysis, it is important to acquire a data set 


42 R.P. Srivastava 


over several orders of magnitude but in a single tectonic environment and rock type 
is necessary. 

Other problems arise if the size of the sampled domain is very small or if the 
domain is so small that it is below resolution limit of the measurement. These 
effects are termed as truncation and censoring after Bonnet et al. (2001). Truncation 
is the underestimation of small fractures due to resolution limitations. These effects 
can be seen in the density distribution. The slope of the power-law curve then goes 
through zero and becomes positive for the smallest fractures. Censoring occurs 
when a limited size of the area is examined compared to the extent of the target 
object. For example, if a large fault being investigated is sampled in a limited area, 
then it can be completely missed or erroneously interpreted. 


8 Some Examples 


Figure 9 is based on the observed data from the faults in one location during a field 
work. It shows the linear relationship between fault zone thickness (fault core) and 
fault throw in a log-log plot. However, if plotted in a log—linear scale it will show 
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Fig. 9 Relation between fault core thickness and fault throw. In a log—log plot this is a linear 
relationship. The data is from outcrop study in a faulted sandstone 
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power-law behavior. Though the observation points are scattered and it is difficult 
to fit a unique line, this kind of observation gives an opportunity to find minimum, 
maximum and most likely cases which are an important input in the fault seal 
modelling. The plot shows that the wider fault zone (fault core) implies larger fault 
throw and vice versa. 

Mathematically, the above observation can be modelled using fractal model of 
the fault throw versus fault thickness relationship given by the following equation: 


log(T) = Nlog(L) —c (2) 


where T is the fault throw; N is the model slope, also known as the scaling 
exponent, which is often a characteristic value in a given region, or location, L is 
fault length and c is the model intercept. Constant c can be determined by fitting a 
line through the observed points. 

The above relationship becomes even more important when there is less data 
available, and we need to rely on some kind of model. To calculate transmissibility 
multipliers of the fault, fault zone thickness (fault core) information is needed, 
which can be computed using above empirical relationship. 

The example shown in Fig. 10 illustrates fault termination and displacement 
profile. Often such instances are visible when you walk on an area with faults, and 
you only see a planer view of the area. In such cases faults can be traced and as 
shown in Fig. 11 a profile of displacement along the fault length can be drawn. This 
information helps to understand the growth history of the faults and helps to 
understand the continuation of the fault when fault through is below seismic 
resolution. 

For instance, in Fig. 10 it may happen that the maximum displacement in the 
centre of the fault is visible in seismic, but as we approach the tip of the fault, the 
throw of the fault gets smaller and smaller, and in most of the cases it is not possible 
to observe the tip of the fault in seismic data. However, during seismic 
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Fig. 10 Model for the fault displacement (slip) distribution. Slip decreases towards the fault tip, 
and maximum in the central part 
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Fig. 11 Displacement (slip) distribution along length profiles in a normal fault 


interpretation, geophysicist may terminate the fault interpretation much before than 
they really exist in nature. This is where such observations in the field help to aid 
the geophysical interpretation of the seismic data, where things are below seismic 
resolution and a geological concept is needed to extend the interpretation. The 
profiles drawn in Fig. 11 give a mathematical estimate of the throw away from the 
maximum observed through, which is a linear relationship. 


9 Conclusions 


Several data- based and experimental studies show that properties of faults follow 
fractal behavior. However, most of the studies on natural faults lack adequate 
collection of data. A better and detailed data collection could help to improve the 
understanding of the fractal behavior of the faults. Given the current status of the 
knowledge available, it is best to use the fractal properties of a fully known fault 
population as an analogue of the properties of an undiscovered fault population. 
Also, with the fractal model of the fault prediction it is important to incorporate 
spatial patterns of available ‘well’ and ‘seismic’ information. The fractal modelling 
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is not a unique modelling, so the results are not deterministic, however, they can be 
converted into a fault probability map. Spatial correlation of fractal based fault 
models can be confirmed from the density of fault-intersecting drill holes mapped in 
the study area. Also, spatial patterns of the fractal based predicted faults can be 
modelled using variograms and checked with the continuity in the spatial patterns 
of known faults. 
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patient even after my never-ending delays in writing this chapter. 
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Detrended Fluctuation Analysis 
of Geophysical Well-Log Data 


D. Subhakar and E. Chandrasekhar 


Abstract Geophysical well-log data provide a unique description of the subsurface 
lithology, as they represent the depositional history of the subsurface formations, 
vis-a-vis the variation of their physical properties as a function of depth. However, a 
correct identification of depths to different lithostratigraphic units is possible only by 
using effective data analysis tools. In the present study, detrended fluctuation 
analysis (DFA) technique has been applied to gamma-ray log, sonic log and neutron 
porosity log of three different wells, A, B and C, located off the west-coast of India 
(i) to discuss the statistical characterization of different subsurface formation prop- 
erties based on their fractal behavior and (ii) to identify the depths to the tops of 
formations by comparing the results of DFA with those of wavelet analysis. The 
DFA technique primarily facilitates to understand the intrinsic self-similarities in 
non-stationary signals like well-logs by determining the scaling exponent in a 
modified least-squares sense. In the present study, DFA was carried out in two ways 
using (1) non-overlapping window method to determine the global scaling exponent 
and (ii) overlapping window method to determine the local scaling exponent. In the 
non-overlapping window method, data segments of different windows, each having 
equal length, were first used to estimate the average fluctuations. The linear 
least-squares regression between the logarithm of average fluctuations and the 
logarithm of window lengths then defines the global scaling exponent. For 
gamma-ray logs of all the three wells, the non-overlapping window method shows 
two distinct ranges of global scaling exponents, in the ranges 0.5—1.0 and 1.0-1.6. 
While the former signifies the presence of persistent long-range power-law corre- 
lations, indicating the stochastic nature of the sedimentation pattern in the data, the 
latter indicates the existence of short-range correlations of non-stochastic nature but 
cease to be of power-law form. On the other hand, the sonic and neutron porosity 
logs of wells A and C show a single global scaling exponent value of greater than 
1.0, signifying the non-stochastic nature of the interval transit time (primary 
porosity) and neutron porosity, respectively, in the entire data sequence as a function 
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of depth. However, in case of well B, the sonic and neutron porosity logs show two 
distinct ranges of global scaling exponents, one in the range 0.5—1.0 and the other 
between 1.0 and 1.5, probably suggesting the effect of different diagenetic conditions 
in well B, compared to those in wells A and C. Choosing a particular window length 
and sliding it with unit shifts over the entire length of data for estimating the 
continuous variation of local scaling exponents as a function of depth defines 
the overlapping window method. This has been applied on all the log data sets of all 
the wells to generate the plots of variation of local scaling exponents as a function of 
depth. Comparison of such plots of variation of local scaling exponents of all the 
logs with the wavelet scalograms of respective logs revealed that the obtained depth 
estimates agree well with the known lithostratigraphy of the study region. 


1 Introduction 


There has been a steady progress in the geophysical well-log data analyses for the 
past three decades, right from the early visual inspection methods to the develop- 
ment of new and fast data analysis techniques to analyse and interpret geophysical 
well-log data for a better understanding of the subsurface formations they offer, in a 
broader perspective (Chandrasekhar and Rao 2012). Principal component analysis 
(Wolff and Pelissier-Combescure 1982), multivariate analysis and non-parametric 
regression (Lee and Datta-Gupta 1999), artificial intelligence (Lim et al. 1999) and 
cluster analysis (Antelo et al. 2001) techniques were used to analyse well-log data 
for detection of electrofacies.' Anxionnaz et al. (1990) carried out cluster analysis 
for the identification of lithofacies. Later, Fourier analysis (Tiwari 1987) and 
semivariogram analysis (Jennings et al. 2000) were also employed to study peri- 
odicities and cyclicities in the data and assess the degree of similarity between 
sample pairs as a function of separation distance within the subsurface formations. 
Kumar and Kishore (2006) studied the classification of electrofacies, by joint 
application of neural networks and cluster analysis. 

Another novel technique, known as wavelet analysis (WA), has been proven to 
be one of the most efficient ones for providing the space localization of different 
subsurface formations. WA has found its wide applications in effectively describing 
the inter-well relationship (Jansen and Kelkar 1997), determining the sedimentary 
cycles (Prokoph and Agterberg 2000), reservoir characterization studies (Panda 
et al. 2000; Vega 2003) and determining the space localization and identifying the 


‘A set of log responses that characterizes a bed and distinguishes it from others (Serra and Abbot 
1980). 
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depths to the tops of formations (Chandrasekhar and Rao 2012). Of late, it has been 
recognized that well-log data display multifractal behavior and thus can be studied 
using fractal analysis to determine scaling exponents, which in turn are related to 
the physical properties of the subsurface formations. Khue et al. (2002) applied a 
generalized multifractal analysis to study the classification of sediment formations. 
Later, Lopez and Aldana (2007) made a wavelet-based fractal analysis of well-log 
data for facies recognition and classified different formations based on their fractal 
dimension. More recently, Hernandez-Martinez et al. (2013) applied detrended 
fluctuation analysis (DFA) to well-log data for identification of facies associations. 
Chandrasekhar and Dimri (2013) described a methodology for wavelet-based 
fractal analysis. 

In the present study, we employ the DFA technique to gamma-ray log, sonic log 
and neutron porosity log data sets of three different wells located off the west-coast 
of India, to (4) statistically characterize different properties of subsurface formations 
based on the fractal behavior of well-log data and (ii) identify the depths to the tops 
of formations by comparing the results of DFA and wavelet analysis. The organi- 
zation of the chapter is as follows. Section 2 briefly discusses the geology of the 
study area. Section 3 details the data used in the present study. Section 4 discusses 
the basic theory of DFA technique and sequential steps to be followed to determine 
the scaling exponents by DFA. Section 5 describes the methodology employed in 
the present work, wherein we have separately used non-overlapping window 
method and overlapping window method to determine the global and local scaling 
exponents, respectively. Section 6 provides the results and discussion of the 
interpretation of the results. Section 7 provides the conclusions of the present study. 


2 Geology of the Study Area 


The Bombay High is the largest hydrocarbon field in India, covering an area of 
about 1500 km?, located in the western offshore basin (Fig. la). The oil and gas 
reservoirs are mostly located in the Cenozoic sediments and are also found in the 
fractured basement, with the most important among the accumulations being those 
in the Miocene carbonate section (Bhandari and Jain 1984). The main oil and gas 
reserves (pay zones) are L-I, L-II, L-II] (limestone) and S-I (sandstone). The 
deepest L-III pay zone (located at about 1300 m) of early Miocene is the largest and 
the most important of all the pay zones, having a thick oil column with a consid- 
erably thick gas cap over the crest. The L-II zone (located at about 990 m) of the 
Middle Miocene is the next important pay zone. Both L-I and S-I are gas reservoirs 
of middle Miocene. Figure 1b provides the geological and lithostratigraphic details 
of the study region. 
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Fig. 1 a Geographical 
location of the study region 
(Bombay High) and its 
contiguous regions in Western 
India (after Chandrasekhar 
and Rao 2012). 

b Lithostratigraphic units of 
Bombay High region (not to 
scale) (after Bhandari and Jain 
1984) 


3 The Database 
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Clay and claystone 
unconformity 
Thick shale with minor limestone 
unconformity 
Limestone with minor shale 
unconformity. 

Shale with prominent fine 
sandstone, siltstone bands 
in the middle 
unconformity: 

Thick limestone with dark grey 
and green shale 
Limestone / shale alteration 
unconformit 
Trap wash 
‘non-conformity 
Basalt / Archean Metamorphics 


For the present study three logs: gamma-ray (GR) log, sonic (DT) log and neutron 
porosity log were used from three wells: well-A, well-B and well-C, located in 
Mumbai offshore (Fig. la), which is one of the largest oil and gas producing basins 
of India. The data were procured from Oil and Natural Gas Corporation (ONGC), 
India. The wells are separated roughly by a distance of about 10 km from each other 
and each represents a vertical section of approximately 500 m below the sea floor. 
The data of all the wells were sampled at an interval of 15 cm. Prior to further 
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Fig. 2 Raw data sets of gamma-ray log, sonic log and neutron porosity log of a Well-A, b Well-B 
and ¢ Well-C of Bombay-High region 
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Fig. 2 (continued) 
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(c) Well C 
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Fig. 2 (continued) 


analysis, the amplitudes of neutron porosity logs of all the wells were magnified by 
a factor of 100 to obtain better statistical inferences. Figure 2 shows different log 
data sets as a function of depth (with depth increasing from top (sea floor) to 
bottom), corresponding to well-A (Fig. 2a), well-B (Fig. 2b) and well-C (Fig. 2c). 
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4 Basic Theory and Estimation of Scaling Exponents 
by Detrended Fluctuation Analysis 


The DFA technique primarily facilitates to study the intrinsic self-similarities in 
non-stationary signals in a modified least-squares sense. By self-similarity we 
mean, the fractal nature of the signal under study. The self-similar nature of a signal 
is mathematically expressed by (Goldberger et al. 2000; Peng et al. 2000) 


x(t) £s’x(t/s) (1) 


where & implies that the statistical properties on both sides of Eq. (1) are equal. In 
other words, a time sequence, x(t), rescaled on one axis, by a factor, s, (t > t/s) and 
on the other axis by a factor of s’(x — s’x), bears identical statistical properties, with 
y being the self-similarity parameter (see Goldberger et al. 2000 for more details). 

Generally, the DFA describes a procedure to determine the fractal behavior in 
the form of scaling exponents of signals generated from stochastic and 
non-stochastic processes. For estimation of scaling exponents through DFA, the 
signal under investigation must be unbounded, which can be converted to a 
self-similar process by integration (Goldberger et al. 2000). The thus generated 
integrated series is divided into short windows (of equal length) and the average 
fluctuation associated with each window of data is obtained by calculating the trend 
(linear least-squares fit) of the selected window and removing it from each data 
point of the corresponding window. This is repeated for various window lengths of 
data. The scaling exponents are then determined from the least-squares linear 
regression between the logarithm of average fluctuation and the logarithm of length 
of the corresponding window. Mathematically, the sequential procedure to be 
adopted in the DFA technique is as follows. 


1. First generate an integrated series y(m) of the N-point spatial data sequence, say, 
x(s), by estimating 


m 


y(m) = S° [x(i) — x] m= 1,2,3...N (2) 


i=l 


x being the average of N data points. As one can easily observe, the integrated 
series thus obtained (Eq. 2) is similar to a random walk and the outcome of each 
value is based on its previous value. 

2. Next, divide the m-length integrated series into various m/k non-overlapping 
windows of equal length, each consisting of k number of samples. 

3. Compute the least-square fit to the data points of each window. This represents 
the local trend y;(m). 

4. Detrend the integrated series, y(m), by subtracting the local trend y,(m) from 
corresponding window. For a window of length k, the average fluctuation, F(k), 
of the detrended signal is calculated by 
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m 


FQ) =|) bd — nO? 3) 


5. Repeat steps 3 and 4 for various window lengths (k) to provide a power-law 
relationship between k and F(k), given by F(k) = k’. Where y denotes the 
scaling exponent, determined by linear least-squares regression between 
log(F(k)) and log(k), signifying the fractal nature of the data. 


Different window sizes (as in step # 2) can be chosen, following the procedure 
described in Peng et al. (1994). Different values of scaling exponents designate 
different physical processes responsible for the generation of the signal under study. 
For example, if a signal is an uncorrelated random noise, then y = 0.5. For 
stochastic processes, varies in the range 0.5—1.0 and for non-stochastic processes 
y > 1.0 (Peng et al. 1994; Hernandez-Martinez et al. 2013). Generally, the scaling 
exponents varying between 0.5 and 1.0 indicate definite power-law correlations 
present in the data. Whereas, the scaling exponent values greater than 1.0 indicate 
the existence of persistent positive correlations in the data but cease to be of 
power-law nature. For more theoretical details and applications of DFA technique, 
the reader is referred to Goldberger et al. (2000). 


5 Methodology 


Figure 3 depicts an example plot of integrated series of gamma-ray log of well-A 
generated using Eq. (2) (see step # 1 in Sect. 4) for a window length of 53 m. The 
linear trend shown in each window depicts the linear least-squares fit corresponding 
to the data segment of that window. To determine the global and local scaling 
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Fig. 3. An example of integrated series (dotted line) generated from gamma-ray log of Well-A 
using Eq. (1). A first-order least-squares fit (solid lines) corresponding to each window of length 
53 m is also shown 
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exponents for well-log data sets of three wells, we have applied two methods, viz., 
(i) non-overlapping window method and (ii) overlapping window method, 
respectively, to the integrated series (Eq. 3). The procedure to estimate the scaling 
exponents by non-overlapping and overlapping window methods is rather different. 
It is described below. 


5.1 Non-overlapping Window Method to Determine Global 
Scaling Exponent 


In this approach, first the data segments of different windows, each having equal 
length, varying in the range 0.75—125 m were considered, following the criteria 
prescribed by Peng et al. (1994) for choosing the window lengths. A first-order 
least-squares fit was fit to the data of each window. The local trend (least-squares 
fit) was removed from each point of the integrated series of the respective windows 
to calculate the average fluctuations, F(k), corresponding to each window (see step 
# 4 in Sect. 4). The estimated average fluctuations of all the windows are averaged 
out to determine the average fluctuation corresponding to that window length. This 
procedure is repeated for various window lengths, till the window length becomes 
equal to N/4, where N is the total number of data points (Peng et al. 1994). Figure 4 
shows such plots of linear least-squares regression fit between log(F(k)) and 
log(k), corresponding to all the well-log data sets of all the wells. The data sets 
show two distinct global scaling exponents, y,; and y,, for gamma-ray logs of all the 
wells and sonic and neutron porosity logs of well-B, indicating short-range and 
long-range autocorrelations in the data respectively. They signify two statistically 
different processes associated with the subsurface formations. Interestingly, the 
boundary depicting the distinction between these two processes is clearly seen in 
the F(k) values between the window lengths, 4.0-4.5 m, for all the wells. The 
presence of correlations in the data could be tested by shuffling the original data and 
comparing the scaling exponents obtained for the original and shuffled data.” If 
y = 0.5 for the shuffled data (i.e. uncorrelated random noise), then the observed 
correlations are due to serial auto-correlations present in the original data. If 
y # 0.5, then the observed correlations are due to the broad probability distribution 
in the data (Kantelhardt et al. 2002). This is tested in the present study. Figure 5a 
shows an example plot of log(F(k)) and log(k) drawn for shuffled gamma-ray log 
data, for which the obtained global scaling exponent, y = 0.5. Figure 5b describes 
the probability distribution of the data, which clearly shows the absence of any data 
beyond the 98 % confidence level of the Gaussian distribution curve. Both these 
tests clearly confirm the occurrence of distinct global scaling exponents and the 
associated correlations are mainly due to the presence of serial auto-correlations 
present in the original data and not due to the broad probability distribution. 


Original data when shuffled, becomes uncorrelated random noise. 
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Fig. 4 Double logarithmic plot depicting the linear correlation between the average fluctuations, 
F(k) (see Eq. 3) and the corresponding window length, k, for all the log data sets of all the wells. 
y, and y, in a, b, c, e and h respectively designate global scaling exponents corresponding to the 
short-range and long-range correlations in the respective log data sets. y in d, f, g and i represent 
the single global scaling exponent of sonic and neutron porosity logs of wells A and C 


5.2 Overlapping Window Method to Determine Local 
Scaling Exponent 


This approach is used to understand the variation of local scaling exponents as a 
function of depth. This helps to identify the depths to the tops of formations in the 
study area, by comparing with those obtained earlier by other techniques. In this 
approach, first a particular window size was considered and the calculated scaling 
exponent corresponding to that window was attributed to the centre of that window. 
Next the window was shifted by a unit step and the above procedure was repeated 
to determine the scaling exponent and attributed it to the centre of that window. 
This is continued till the end of the data is reached. We call the scaling exponent 
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Fig. 5 a Global scaling exponent (y) value determined for the shuffled gamma-ray log data of 
well-C. The shuffled sequence has resulted in an uncorrelated random noise, whose scaling 
exponent value is 0.5. b Probability density function and the Gaussian normal fit to the gamma-ray 
log of well-C. Note the absence of fat tail distribution in the data 
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Table 1 Comparison of the depths to the tops of formations delineated from the plots of variation 
of local scaling exponents as a function of depth (Figs. 6, 7 and 8) with those obtained from 
wavelet analysis using Gaussl wavelet for different well-logs (Chandrasekhar and Rao 2012) 
Known depth estimates provided by ONGC Ltd. are also given. Note the good agreement between 
the observed and estimated values 


Reservoir top ONGC Lid. Wavelet DFA (present study) 
analysis 
(Gauss 1) 


GR 


thus obtained as a local scaling exponent. The entire procedure described above is 
repeated for different window lengths, varying in the range 16-36 m 
(Hernandez-Martinez et al. 2013), which corresponds to 4 to 8 times the cross-over 
window length of 4.0-4.5 m (cf. Sect. 5.1), depicting the boundary between 
short-range and long-range correlations in the signal. Figures 6, 7 and 8 show the 
plots of variation of local scaling exponents as a function of depth determined for 
all the logs of wells A, B and C, respectively, corresponding to a window length of 
23 m. This window length was found to be optimum for clear identification of 
depths to different subsurface litho units, when compared with those of the wavelet 
analysis results of same data sets (Chandrasekhar and Rao 2012). Table | provides 
a comparative description of the depths to the tops of formations obtained from 
wavelet analysis, from ONGC Ltd. and the present study. 


6 Results and Discussion 


Figure 4 shows the global scaling exponents derived for gamma-ray log, sonic log 
and neutron porosity log of all the wells using non-overlapping windows of various 
sizes. The two global scaling exponents, y,; and y,, shown in Fig. 4 correspond to 
two different groups of window lengths varying in the ranges 0.75—-4 m and 
4.5-125 m, indicating short-range and long-range correlations, respectively, in the 
signal. Gamma-ray logs of all the wells (Fig. 4a—c) and sonic and neutron porosity 
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Fig. 6 Variation of local scaling exponent as a function of depth, calculated for all the data sets of 
Well-A using an optimum window length of 23 m. Note the clear demarcation of depths to the tops 
of L-I, L-I, L-III and S-I zones, identified by the DFA technique 


logs of well-B (Fig. 4e, h) show two distinct global scaling exponents, y, and yp, 
being greater than | and less than 1, respectively, indicating the presence of 
non-stochastic and stochastic processes in the subsurface formations (cf. Sect. 4). 
However, the sonic and neutron porosity logs of wells A and C (Fig. 4d, f, g, 1) 
show a single global scaling exponent value of greater than 1, signifying the 
non-stochastic nature of the interval transit time (primary porosity) and porosity, 
respectively, in the entire data sequence. 

As is well known, gamma-ray logs facilitate to understand the subsurface sedi- 
mentation patterns due to the large differences in the natural radioactivity of shaly 
and non-shaly formations. Thus the two distinct global scaling exponents deter- 
mined for gamma-ray logs of all the wells (Fig. 4a—c) clearly demonstrate the distinct 
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Fig. 7 Same as Fig. 6, but for Well-B 


fractal behavior of sedimentation pattern, such that it appears to be non-stochastic 
when window lengths are small and stochastic when window lengths are large. 
Sonic logs, which record the interval transit times, are mainly used to estimate 
the primary porosity of the formations and neutron porosity logs, which measure 
the hydrogen index of the formations, provide true porosity in case of limestones 
(Serra 1984). The single global scaling exponent values of greater than 1.0 derived 
for sonic and neutron porosity logs of wells A and C (Fig. 4d, f, g, i) indicate the 
nature of the primary porosity and the true porosity in the entire formation sequence 
to be non-stochastic. It is interesting to observe that the fractal behavior of both 
these logs does not show any significant differences in the global scaling exponents 
although they are obtained with different window sizes (see Fig. 4d, f, g, i). This 
suggests that the differences in the successive F(k) values are very small and thus 
there is a steady increase in them with increase in window length. Usually, this 
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Fig. 8 Same as Fig. 6, but for Well-C 


happens when there are no significant changes in the petrophysical properties of the 
formations within the considered depth range. Thus, the estimated global scaling 
exponents for both the logs of wells A and C confirm the non-stochastic trends in 
the primary as well as true porosity. Since the primary porosity generally decreases 
with depth due to compaction and diagenesis (Medina 2009), its non-stochastic 
trends seen in the data corroborate this observation. 

In case of well-B, both sonic and neutron porosity logs, unlike in the cases of 
wells A and C, show two distinct trends with respect to window lengths, resulting in 
two distinct global scaling exponents (Fig. 4e, h). Also, in case of gamma-ray log of 
Well-B, the difference between y, and y, values of Well-B (Fig. 4b) is slightly 
higher than those observed for the same logs of well-A and well-C. Thus the overall 
subsurface nature of well-B appears to be quite different from that of wells A and C. 
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Perhaps, this could be due to the effect of different diagenetic conditions in well-B. 
However, more log data sets need to be studied from this area to further confirm this 
observation. 

Figure 4a—c, e, h, show a change in global scaling exponents in the cross-over 
window length of 4.0—4.5 m, signifying the boundary between non-stochastic and 
stochastic behavior in the data of respective logs of all the wells. 
Hernandez-Martinez et al. (2013) suggest that for calculating the local scaling 
exponents as a function of depth using overlapping window method, the window 
length should be four times that of the cross-over window length. They opine that 
the local scaling exponents thus calculated, facilitate to demarcate the subsurface 
lithounits. Accordingly, various window lengths starting from 16 to 36° m were 
chosen to calculate the local scaling exponents as a function of depth for all logs of 
all the wells. A comparison of all such plots generated using various window sizes 
with the wavelet analysis results of well-log data sets of the same region (see 
Figs. 5-7 and Table | of Chandrasekhar and Rao 2012) clearly showed that the 
window length of 23 m is optimum to identify the depths to the tops of reservoir 
zones in all the wells (Figs. 6, 7 and 8). Beyond the window length of 23 m, the 
local scaling exponents are largely averaged out and did not provide proper reso- 
lution to identify the depths to the tops of the reservoir zones. 

A careful examination of the plots of local scaling exponents obtained as a 
function of depth for all logs of all the wells (Figs. 6, 7 and 8) and Table | clearly 
identifies not only the depths to the tops of L-I, L-II, S-I and L-III reservoir zones, 
but also the probable locations of unconfirmities between L-II and L-III zones (see 
Fig. 1b). A careful observation of Fig. 6 shows the sudden changes in local scaling 
exponent values in the range 0.8—1 (between 1090 and 1320 m) to I-1.3 (1320 m 
and below). These boundaries at 1090 and 1320 m clearly depict the plausible depth 
location of unconfirmities (see also Fig. 1b). The likely depth locations of uncon- 
firmities in wells B and C are also shown in Figs. 7 and 8. 

Limestones show low gamma-ray counts due to the absence of radioactivity. 
However, since L-I, L-I and L-III zones depict the limestone-shale boundaries, the 
variation of local scaling exponents of gamma-ray log with depth corresponding to 
these boundaries indicates positive high (Figs. 6, 7 and 8), due to increase in the F(x) 
values as the window approaches the shale boundary. Similarly, in case of sonic log 
since the interval transit time is low in limestone reservoirs and is high in shales, the 
limestone-shale boundaries marked by positive high local scaling exponents could be 
easily identified (Figs. 6, 7 and 8). In case of neutron porosity log, the hydrogen index 
is high in shales due to the presence of bound water compared to hydrocarbon-bearing 
limestones. Accordingly, the local scaling exponents show positive high value at the 
limestone-shale boundaries as a function of depth (Figs. 6, 7 and 8). 

In case of S-I (the boundary depicting the sand and shale) gas-bearing zone, 
since the gamma-ray log shows gradational changes due to mixture of sand and 
shale, the identification of such boundary became difficult by the DFA method, as 


This maximum window length corresponds to 8 times the cross-over window length of 4.5 m. 
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has been the case for gamma-ray log of all the wells. In gas-bearing zones while the 
interval transit time as seen in sonic logs will in general be higher compared to that 
of non-hydrocarbon zones, the hydrogen index as measured in neutron porosity 
logs is generally low compared to that of non-hydrocarbon zones. Accordingly, the 
S-I zone could be clearly identified in wells A and B (Figs. 6 and 7). However, the 
S-I zone could not be clearly delineated in well-C (Fig. 8) and the reasons for this 
could not be ascertained in the present study. 


7 Conclusions 


DFA is proven to be one of the efficient mathematical techniques for studying the 
fractal nature of non-stationary signals like well-logs. Among the non-overlapping 
window method and overlapping window method employed to determine the global 
and local scaling exponents, the former has proven to be an efficient method to 
clearly distinguish the non-stochastic and stochastic trends in the physical proper- 
ties of the formations and the latter to identify the depths to the tops of reservoir 
zones. Particularly, the distinct stochastic and non-stochastic behaviors observed in 
the analysed data sets could confirm the nature of the subsurface formations in 
well-B to be quite different from those of wells A and C. However, we feel that a 
study on more data sets from this region might strongly support our observation. 
For the first time, through this study, we have made a comparative study of wavelet 
analysis results with those of the present study and identified a window length of 
about 23 m to be optimum to demarcate the depths to the tops of reservoir zones 
and showed that the results obtained from the DFA technique and wavelet analysis 
match quite well. Further studies using wavelet-based fractal analysis and multi- 
fractal DFA, applied to various well-log data sets of diverse regions, would be 
worthwhile for a comprehensive understanding of the nature of the subsurface 
formation properties in a broader perspective. 
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Fractal Characters of Porous Media 
and Flow Analysis 


Pallavi Banerjee Chattopadhyay and Nimisha Vedanti 


Abstract Porosity is a complex multivariate function which controls fluid flow 
through porous media. The variables essentially are the characteristics of pore 
structures such as type, size, shape and arrangement of pores; pore space connections; 
area of pores that is open for flow; tortuosity of the flow paths and composition of 
pores, etc. Predicting the rate of flow and the flow patterns of fluid in bulk geological 
formation is very important for many environmental problems and oil industry. 
Fracture networks are secondary porosity that is known to exist within the subsurface 
geology that is expected to influence the flow through geological heterogeneous 
irregular porous media. Many of them are relevant to the migration and entrapment of 
fluids within the reservoir. It shows the nature of symmetry in geological formation. 
However, no general framework exists to systematically study the fluid flow through 
the fractured subsurface geometry which is complex in nature. Since direct mea- 
surements of flow through complex permeable media are time consuming and require 
experiments that are not always feasible, an analytical model could be more useful. 
Therefore, the focus of this chapter is on the development of simple analytical models 
based on medium structural characteristics to explain the flow in natural fractures. 
The pattern of fracture heterogeneity in reservoir scale of natural geological forma- 
tions can be viewed as spatially distributed self-similar tree structures. Application of 
fractal geometry is useful to define the porous structure, where fractal geometry is the 
study of mathematical shapes which display self-similar, meandering and tortuous 
porous media details. Based on the above understanding, a fractal model of con- 
tinuum percolation is presented here, which quantitatively reproduces the flow 
path geometry by applying Poiseuille’s equation, where flow in fractures is driven by 
the pressure differences at the two ends of the path. 
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1 Introduction 


Interaction between fluids and pore structures delineates the fluid motion in porous 
media. Therefore, structure of fractures within the geological formations and its role 
in fluid flow is an important topic of research. Naturally fractured porous media are 
frequently encountered in subsurface geological formations that include ground- 
water, hydrocarbon and geothermal reservoirs (Moench 1984; Chang 1993). Matrix 
(primary porous formation) usually has low permeability and high storability, while 
fractures have high permeability but low storability. This implies that flow through 
the formations is highly dependent on the distribution pattern of fractures and their 
interaction. Many researchers (Barenblatt et al. 1960; Warren and Root 1963) have 
used double porosity concept in reservoir simulation where two different partial 
differential equations define matrix and fractures flow. Various researchers (Ueda 
et al. 1989; Lim and Aziz 1995; Sarma and Aziz 2006) have attempted to improve 
the basic flow equation by using various approaches that account for fracture flow 
with practical importance. Predicting flow through a fractured rock is one of the 
most challenging and difficult problems. In large systems such as reservoirs, it is 
challenging to obtain sufficient geometric and hydraulic information which con- 
strain the connectivity of fracture networks and define fluid flow trajectories. 

Subsurface fracture distribution is a non-linear system in nature where magni- 
tude and direction vary in the wildest way. During the 1980s, physicists and some 
hydrogeologists had put in considerable amount of effort to understand the physical 
properties of fluid flow in fractured porous media (Sen et al. 1981; Krohn and 
Thompson 1986). Petroleum engineers and hydrologists generally deal with con- 
solidated fractured rocks where the conceptual structure of the medium is critical 
for estimation. Probably the most striking aspect of naturally fractured porous 
media is the complexity in texture. Fracture sets typically occur as groups of tens to 
thousands of individual fractures in reservoirs and obtaining a useful and valid 
description of the media is the fundamental difficulty in estimating flow and 
transport properties. It is a common practice to narrow down the complexity by 
largely ignoring the geometrical complications and many researchers (Sudicky and 
McLaren 1992; Bogdanov et al. 2003; Helmke et al. 2005) have represented such 
matrix as a discrete fractured network to focus on the core problem at hand. The 
most important reason for modelling fractures discretely is to represent the con- 
nectivity of flow paths in the subsurface using a geologically realistic and simple 
heterogeneous representation. The study of discrete fracture networks is more 
important in systems where fractures are large and sparse compared to the scale of 
reservoir and in locations where fracture geometry can be characterized. However, 
such characterization may be limited to the large spatial extent of a reservoir. 
Commonly fractures are developed in consolidated rocks, where the fractured 
medium forms a preferential flow path over the primary porous path. Fractures 
control fluid flow in the geological settings which are anisotropic and heteroge- 
neous systems, characterized by a network of fractures spread over various regions 
of porous matrix in a large-scale reservoir. 
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Katz and Thompson (1985) noticed that most of the porous media are fractals, 
where steady-state crystal growth during rock formation is a plausible cause for the 
self-similar geometry. Jacquin and Adler (1987) performed an extensive series of 
numerical simulations and the results show that the permeability of the porous 
media can be approximated by a power-law function of the porosity. Prediction of 
fluid flow through complex fractured systems as fractures spread over a wide range 
of scales is very complicated (Barton et al. 1983; Brown 1989) and a traditional 
approach cannot demonstrate quantitative characteristics at reservoir scale. Many 
conventional models of fracture networks have been developed as a planar system, 
whereas natural fractures are not generally planar (Gertsch 1995). Many researchers 
(Mandelbrot 1982) have represented fractures as a finite line segment. However, the 
effect of internal surface roughness of the fractures is not considered in most of the 
fracture flow models. Snow (1969) proposed the simplest model of planar fracture 
that is the space between two parallel and infinite flat planes. First, Wang et al. 
(1988) used fractal geometry to overcome this difficulty. He proposed a fractal 
model of fractures which considered the aperture corrections. Reservoir scale 
fracture distributions commonly exhibit similar features at different resolutions 
(Fig. 1). 

Most of the porous media models are based on single pores, however averages 
are performed over a pore-size distribution without considering the effects of 
connectivity (Kozeny 1927; Carman 1956; Mualem 1976). More complicated 
averaging schemes, such as Burdine’s (1953), which uses a joint probability dis- 
tribution rather than a single probability, diminish the connectivity of the large and 
wide pores in an area. 

Thus, interest of the present study is to explore the fluid flow behavior within a 
fractured geo-fluid reservoir. Multiscale fault and fracture networks represent the 
most relevant features allowing geo-fluid migration within reservoir scale subsur- 
face structures. In fractured reservoirs, the largest faults are usually detected by 
seismic survey and smaller faults are captured by well logging. However, field 
survey can only provide very limited information about fracture clustering and 
connectivity, which are all together responsible for flow properties. Therefore, 
reservoir scale study demands realistic flow models. 


Fig. 1 Fractal geometry can relate the fracture distribution at one scale to other (Martin et al. 
1998) 
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Dimri et al. (2012) provide the fundamentals of fractal applications in reservoir 
geophysics and reservoir simulation study. Reservoir scale fracture distribution is a 
self-organized process on geological time scale (Bak et al. 1987; Cowie et al. 1993). 
Natural fracture networks show a power-law distribution (Walsh and Watterson 
1992; Cowie et al. 1995). Natural display of fracture distribution within a reservoir 
shows fractal-like property. In the general sense, fractal objects show the properties 
of being nearly the same at every progressive scale. However, in mathematical 
definition, no natural object is purely a fractal, instead it can be called as an 
‘approximate fractal’ or ‘statistical fractal’ that displays ‘self-similarity’ and 
‘self-affinity’ over extended but finite scale of ranges (Bovill 1996). Fracture dis- 
tribution in reservoir scale is a fractal geometry that is self-similar in pattern and 
highly irregular that can be compared with trees with their branches (Fig. 2). 

While proper attention is drawn towards the fractal geometry and fractal 
dimension, analysis of flow through the geometry is neglected. Fractal model, 
substituting a complex fractured porous media with a network of interconnected 
pores, is a powerful numerical method for flow modelling in porous materials. The 
model is accepted because fractal geometry is developed to describe irregular 
natural shapes and to quantify the scaling properties of fracture networks. 


Fig. 2 Fluid flow through IN FLOW 


saturated pore throats 


OUTFLOW 
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2 Fluid Flow Through Fractal Geometry 


The fluid flow equations for fractured system are very complicated as the flow 
depends on viscosity of the fluid and the diameter of the fractured path. The 
assumption of constant fluid viscosity in the fractured geometry is necessary to 
simplify the flow as a near Newtonian fluid. However, walls in a naturally fractured 
system are irregular and can change their diameter, which can affect fluid flow in 
the media. It is assumed that the fractures have constant diameter for small segment, 
which is necessary for application of hydrodynamic equations and analytical cal- 
culation of the model. 

Turbulent flow is possible in large diameter fractures. When the fluid moves 
slowly, without much mixing between the layers the flow estimation then assumes 
laminar flow for the entire system. Assumption of laminar flow is valid for small 
diameter fractures, distributed in a large area. Such flows are mainly quantified by 
Darcy’s formula. However, Darcy’s equation has serious limitations for fractured 
laminar flow, but it gives the correct result for creeping filter flow through relatively 
long, uniform and isotropic porous media of sand or unconsolidated formations. 
Therefore, Darcy’s law is not suitable for fractured media flow, which has com- 
plicated topology embedded in three dimensions. Three dimensions of fractured 
media are obtained through stereo projection, which may not be always simple. 

Similar to Darcy’s law, Poiseuille’s law is based on the experiment of flow, but 
flows are through small diameter pipes. The proposed flow equations by Darcy 
(Hagen 1839) and (Poiseuille 1840) are based on geometric differences between a 
sand column and small diameter pipes, respectively, for laminar flow. The fracture 
flow is assumed to be laminar through entire capillary size fractal geometry. The 
most important feature of a fractal object is its “degree of irregularity”. The degree 
of irregularity (D) in fractal geometry allows us to quantify the roughness and 
irregularity of a shape in numerical value. Researchers (e.g. Katz and Thompson 
1985; Turcotte 1970; Bayles et al. 1989) suggest that, in a given length scale range, 
natural porous media can be described as a fractal of dimension D. 


2.1 Porosity and Fractal Media 


The fractal fragmentation model of Turcotte (1986) added legitimacy to the study of 
fractal models of porous media, particularly soils, because it developed a mecha- 
nism by which scale-independent fracture properties could form a fractal distri- 
bution of particles. As an extension, Rieu and Sposito (1991) then developed a 
model, named as the Rieu and Sposito (RS) model, based on fractal pore space to 
predict pressure—saturation curves from particle-size data. In the model, d,, and ry, 
denote the maximum pore diameter and radius, and dp and rp the smallest diameter 
and radius. The reason for their choice is that they can use an index i, representing 
the iteration from ‘m’ number to 0 of the fractal process. 
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For simplicity, V; represents the total volume in all pore sizes greater than d,, and 
less than or equal to d;. There is a constant ratio N of the number of pores of 
diameter d;,; = gd; to the number of pores of diameter d;. Here, g is the ratio of pore 
diameters in successive classes and is less than 1. Define the partial volume 
P; = V; — Vis, the total volume as Vo and the volume of the solid material as V,,. 


m 


Vo= ¥ > NVi41 +P; (2.1.1.1) 
i=0 


Here, pores smaller than d,, can be ignored; that is, if we could resolve smaller 
pores, the iteration would proceed further. The total pore volume, V,, can be written 
as: 


1 — (Nq?)” 


Vp =P 
pT Ne 


Oi) 


One can express the porosity and substituting the Pp value (Py = Vo (1 — Nq’)) to 
find: 


& =1-—(Nq*)” (213) 


The numerical factor q,,, however, is nothing more than the ratio of the smallest 
pore diameter to the largest, d,/do, so that fractal dimensionality: 


“log __log(l- 4) _, (t=) ee 
=— = = ALA 
log; log () do 


The above equation relation establishes the RS result for the porosity. The 
power-law distribution of pore sizes is bound by a maximum radius r,, and trun- 
cated at the minimum radius 7p. Knowledge of ¢, ro and 7, is sufficient to explain 
dimension (D) explicitly. 


2.2 Poiseuille’s Equation: Simplified Navier-Stokes 
Equations 


The scaling behaviors are generally attributed to medium heterogeneity 
(Schulze-Makuch and Cherkauer 1997; Carrera 1993). The parameters for fracture 
geometry are functions of the scale, where different scales could be considered 
according to the specific problem (i.e. commonly oil fields are regional scale 
problems). The regional scale fracture distribution supports self-similar assumption, 
very similar to the tree shape fractal model (Rovey 1994). In fractured framework, 
the porosity, which could be considered as the small diameter tube with tortuosity, 
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plays a fundamental role with regard to the water flow in the porous medium. Flow 
occurs from locations with high potential energy to locations of lower potential 
energy in pursuit of equilibrium state. The driving force for flow is called potential 
gradient, the difference in potentials between two points in a system separated by a 
certain distance ‘L’. In Eq. 1, i is the potential gradient and Y; and > are potential 
energies at inlet and outlet, respectively. 


f= (7) (2.1.2.1) 


The gradient is developed due to differences in pressure, position in gravity field, 
chemical concentration, temperature and position in electrical field, which leads to 
spontaneous flow of mass or energy. 

A Newtonian viscous fluid of constant density is shown to be in steady motion 
down a cylindrical tube of radius of R (Fig. 2). The pressure gradient should be 
down the tube and velocity profile is parabolic, because of viscous stresses at the 
walls of the tube. The flow is considered to be steady and laminar of Newtonian 
fluid with constant viscosity through a small diameter of cylindrical tube similar to 
fractured segment. For such a model Navier-Stokes equations are simplified to 
Poiseuille’s equation: 

The Navier-Stokes equation in the coordinate is given as: 


- ue Vue 4 - =u, (2.1.2.2) 

os bu Vay Mt 4 LP = ou, > an) (21.23) 
aug ue Duy te 1” _ »( Ku - i) (2.1.2.4) 
= 1 on cp sau _¢ (2.1.2.5) 


Here, in Eqs. (2.1.2.3) and (2.1.2.4), K= V7 - + and u-V= 
uz a + up & + we It is assumed that the axis of the tube is the z-axis, r is the 
radial variable and u = (u,, u,, Ue) = (u(r), 0, 0) is the velocity field in cylindrical 
polar coordinates. 

Consider a steady laminar flow of Newtonian fluid with constant viscosity 
through a horizontal cylindrical rigid tube (fracture segment). For such a model, 
Navier-Stokes equations are simplified to Poiseuille’s equation. 
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(2.1.2.6) 


Poiseuille’s law relates to the fluid flow Q [ml/s] through a fracture segment with 
the difference in fluid pressure Ap at the two ends of segment. Segment radius is 
represented by 7, length is L and viscosity is y of the fluid. 

The most effective factor controlling flow is radius of the segment. High pressure 
can be caused by narrowing segment and is reduced by increasing the fracture 
radius. During laminar fluid flow, cylindrical liquid in segment is exposed to 
internal friction, which represents the resistance of flow R. 


— 8ub 


a 


(2.1.2.7) 


The resistance of fluid motion through fractured segment is most strongly 
dependent on radius, with the fourth power relationship. Because of fluid friction, 
the velocity of fluid within the fractured segment varies from none in wall proximity 
to maximum value in the centre of the vessel creating a parabolic velocity profile. 
Average velocity (V,) with respect to cross-section is determined by (Fig. 3) 


— Apxr? 
~ 8uXxL 


(2.1.2.8) 
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Fig. 3 Flow rates, derived by power law for geological fractured media, are plotted against 
average diameter of the fractures segment 
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3 Concluding Remarks 


It is essential to be able to predict the rates of flow and the flow patterns of fluids in 
bulk rock during recovery of oil, or for the isolation of nuclear or toxic wastes. 
However, it is challenging to define average flow properties of fractured networks. 
In this chapter, multiscale nature of fracture networks distribution and the conse- 
quences of flow property are discussed. The fracture length distribution could 
reasonably be modelled by a power law whose exponent quantifies the scaling 
nature of fracture networks. The fact that fractures exist at all scales, from micro- 
cracks to pluri-kilometric faults, is a challenging issue for defining average flow 
properties of fracture networks. The fracture distribution in reservoir scale is rea- 
sonably modelled by a power law whose exponent quantifies the scaling nature of 
fracture networks. It has been measured that three-dimensional fracture outcrops 
with values ranging from 2 to 3. The transmissivity distribution is much lesser 
known because of the difficulty in measuring transmissivity in complex hydroge- 
ological sites. Log-normal functions or power laws are classically used as trans- 
missivity distribution models, but this has to be taken as a conjecture. The flow 
model is highly dependent on the scaling property of fracture networks and more 
specifically on the power-law length exponent. The two end-member models are the 
percolation-like model, where all fractures have a length much smaller than the 
system size. The probability of occurrence of large fractures (i.e. larger than the 
studied system) increases with scale, as such systems are always connected above a 
critical connecting scale that is ignored in this chapter. Permeability, as well as 
other flow properties, is closely related to the flow structure, which can be char- 
acterized by length scales, which depend on the fracture distribution and on the 
transmissivity distribution per fracture. Here, in the proposed models of fractal, 
tree-shaped distribution of fractures is accepted for large reservoir scale distribution. 
It has been scaled in such a way that the principal statistical characteristics are 
quantitatively comparable with measurements of real fractures. Possible simplifi- 
cation and generalizations assumed for a flow model for defined asymmetrical 
tree-like distribution have been discussed. Proposed approach enables straightfor- 
ward analytical evaluation of flow parameters, which usually involves very 
sophisticated numerical computational methods especially for such complicated 
phenomena and fracture flow in general. 
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unto Porous Solids 
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Abstract The conventional approach for fractal dimension (FD) estimation using 
box count method has been widely used in the analysis of imageries especially in 
the domain of earth systems modelling and has been known to provide insight into 
the complexities of the system as well as the dynamics of the processes involved. 
However, for heterogeneous imageries such as micrographs, etc., the information 
provided by estimated FD seems to be limited. The present work establishes this 
limitation in the use of FD (using HarFA 5.5 software) and extends the concept of 
fractal dimensioning into lower scale segregation levels and evaluating their dif- 
ferential scores. In this approach, fractal differential adjacency segregation (F-DAS) 
scores are estimated using MATLAB 14.0 for each of the image pixels (of SEM 
imageries) using the arithmetic means of the grey levels of the adjacency pixels 
enclosed by the box (used for counting in the conventional methods). The present 
analysis provides a better understanding of the variability of the system (in this 
case, adsorbents—adsorbate interactions), unexplored by qualitative analysis of SEM 
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imageries as well as the functional groups using FTIR. This work provides sys- 
tematic steps of estimation of F-DAS scores of any imagery, the assumptions 
underlying the approach as well as the scopes of its applications in analysis of 
various earth systems. 


1 Introduction 


Scanning electron microscope (SEM) has been widely used for characterization of 
solid materials at micro-scale and hence has been widely used to detect and analyse 
surface fractures, provide information on microstructures, examine surface con- 
taminations, reveal spatial variations in chemical compositions, provide qualitative 
chemical analyses and identify crystalline structures and so forth. SEM is a 
microscope that uses electrons rather than light to form an image. There are many 
advantages to using the SEM instead of an ordinary microscope. The SEM has a 
large depth of field, which allows a large amount of the sample to be in focus at one 
time and produces an image that is a good representation of the three-dimensional 
sample. The SEM also produces images of high resolution, which means that close 
features can be examined at a high magnification. 

The combination of higher magnification, larger depth of field, greater resolution 
and compositional and crystallographic information makes the SEM one of the 
most heavily used instruments in research areas and industries, especially in the 
semiconductor industry. Since the development of SEM (by Dr. Charles Oatlev in 
the 1950s), its area of application has been increased extensively, especially 
because of its ability for detailed three-dimensional and topographical imaging, ease 
of operation, fast analysis and availability of imageries in digital form (Joseph 
2003). 

However, one of the major disadvantages associated with SEM-analysis is that it 
provides more or less a qualitative understanding of the system, rather than a 
quantitative estimation. Hence, the interpretation of SEM imageries is mostly left 
unto the skill and experience of the interpreter. In fact, there has been a growing 
need for the objective and quantitative estimation assessment of SEM imageries, 
especially with its increasing diversity of application (Price 2005). 


1.1 Fractal Dimensions from SEM Imageries 


Of late, fractal conception of matter developed and explained by B.B. Mandelbrot 
in his book “Fractals Geometry of Nature” (Mandelbrot 1983) has been increas- 
ingly used to understand the self-similar and scale-independent geometry of natural 
structures, where fractal dimension (FD) strictly exceeds topographical dimensions. 
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Because of partial correlation existing among a wide range of natural systems, FDs 
have been used to define a wide range of biological and geological systems. 

The concept of ‘self-similarity of surface at different scales’ is mostly used in the 
fractal analysis of earth sciences data (Dimri 1992). Basically the structure is 
characterized by single descriptor, the FD ‘D’. The FD lies within the range 
2 < D <3, where D = 2 represents a smooth surface and an increasing value of 
D represents an increasing surface roughness/porosity (Othman et al. 2006; Maus 
and Dimri 1994, 1995, 1996). The concept has been used by Dimri (2000) for 
studying the flow media. In fact, FD can also be correlated with various surface 
roughness parameters (Chauvy et al. 1998; Fang et al. 2003; Seitavuopio 2006 and 
Dimri 2000). Bigerelle and Iost (2003) have revealed in their studies that statistical 
comparison of relevant roughness parameters of FD is the most relevant parameter 
to describe the surface topography. The FD of a certain surface can be in principle 
inferred in different ways. These include, but are not limited to: scattering, 
adsorption of gases, impedance spectroscopy and surface image analysis. The latter 
two methods are especially widespread and used in different areas of research 
(Lasia et al. 1999 and Li et al. 2003). 

Among the many methods available for analysis and description of surface 
topographies, the SEM and the atomic force microscopy (AFM) are widely used for 
surface imaging and characterization. Due to its high depth of focus, the SEM can 
provide detailed topographical information about the surface, but cannot provide 
quantitative topographical information. Such information can be obtained by fractal 
analysis of the images revealing the characteristic FDs (Risovic et al. 2008). 

Since the surface dynamics of a solid object, caused by physico-chemical 
stresses, can be captured by the variability of their surface topography in SEM 
imageries, and such effects are largely guided by the surface characteristics of the 
object concerned (porosity, texture, micro-fractures, etc.), which are mostly natural, 
it is legitimate to explore fractal signature of surface dynamics of different mate- 
rials, effected by the subjected stresses (such as sorption, reaction, activation, and so 
forth). In aqueous adsorbents (the solid/liquid that adsorbs the adsorbates, i.e., the 
liquid/gases on their surfaces) is characteristic of the microstructure of the adsor- 
bent, which is non-homogenous and self-similar, FD may be a useful tool for 
characterization of the SEM imageries of the adsorbents. There have been few 
reports with regard to fractal applications in adsorption systems, however, the 
present work is taken up to access the scope of developing a characterization tool 
for various adsorbents, using SEM imageries (Brouers and Sotolongo 2006). 

Of the various methods used in estimating the FD of an imagery, one of the most 
common methods is box-counting method, wherein the grey levels of the image 
have been computed, based on BW (black to white) boundary, B—BW (Black to 
BW) boundary and W—BW (White to BW) boundary. These analyses throw picture 
unto the fractal nature of associated phases (bright/dark) as well as their blends. The 
present work also carries out these analyses to estimate fractal characteristics of the 
imageries, using random samples to account for wide textural inhomogeneity in the 
original imageries itself (details in experimental). However, the FD of the images 
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does not capture the regional variation in complexity, for studying which a novel 
approach has been proposed herewith, namely—fractal differential adjacency seg- 
regation (F-DAS) Scores. 


2 Fractal Differential Adjacency Segregation (F-DAS) 
Scores—Concept and the Approach 


F-DAS scores derive their basis on the estimation of FD, similar to box-counting 
method, but here, instead of counting the boxes, the average of the scores of the 
boxes is used to replace the box values. Thus, with increase in the box sizes, the 
scores of the boxes show changes, based on the relationship of the adjacent scores. 
For example, if a white pixel is surrounded by all black pixels, the average grey 
level of the box containing this white pixel and its adjacent pixels would be close to 
the score of the black pixels. In fact, the change in the value of the particular pixel, 
for any fixed box size, would depend on the grey levels of the pixels constituting 
the box, under consideration. Then, with increasing box size, the value would 
further be affected by the grey levels of newly added adjacent pixels. Thus, the 
underlying principle of box-counting strategy of conventional FD estimation would 
be replaced by a harmonic aggregation (to account for adjacency effect), instead of 
the arithmetic summation. 

The standard power law for FD estimation for box-counting method, i.e. the 
negative of the slope between logarithm of the estimated area versus logarithm of the 
box dimensions used to estimate the area, was also used in this study where the slope 
was calculated for each pixel with regard to its change in values over the change in box 
sizes. Hence, the estimated negatives of the slopes form a matrix, instead of a single 
number, as in the case of conventional FD, thus providing a spatial understanding of 
the complexities associated with the image (on account of associated adjacency). 
Now, for an understanding of the change in the slope, the difference between slopes 
with first two boxes and that between second and third boxes was computed, thereby 
giving scope for enhanced understanding of the adjacency effect with regard to 
persistency of fractal estimates. The resultant matrix is named as F-DAS scores, which 
can be used to evaluate the estimate of persistency against the adjacency effect. 

However, it is important to note that, since this approach deals with moderation 
of the scores based on the adjacent scores, which is non-cumulative (unlike the 
approach used in conventional FD estimation), these values do not obey the 
characteristics of FDs (i.e., non-negativity and limitation to maximum allowable 
dimension, i.e., less than 3 in this case, as expected to be in case of two-dimensional 
features like images). 

In the present study, the F-DAS scores matrices were estimated using the grey 
levels of the various SEM-images and used for understanding the characteristics of 
the adsorbents (i.e., the porous solids acting as the base for adsorption in aqueous 
media). 
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3 Fundamentals of Adsorption and Relevance of Fractal 
Studies 


Adsorption is the process of concentration of one or more contaminants (of either a 
gas or a liquid) on the surface of a micro-porous solid. The solid substance, which 
attracts pollutants, is called Adsorbent, and the pollutant concentrated at the surface 
of the solid is Adsorbate. Adsorption phenomena can be physical, chemical and/or 
biological in nature. Physical adsorption (also called, physisorption) is mainly due 
to van der waals forces and electrostatic forces between adsorbate molecules and 
atoms of adsorbent and is directly proportional to the amount of available solid 
surface. It takes place in one or more molecular layers, often accompanied by 
capillary condensation within the pores (Srihari and Das 2004). This process is 
relatively rapid and reversible. In fact, by lowering pressure and raising the tem- 
perature the adsorbed pollutant can be desorbed (or recovered) in many cases. The 
chemical adsorption (also called, chemisorption), on the other hand, involves a 
chemical bond between the adsorbent and the pollutant molecules (or adsorbate), 
held strongly to the surface of the adsorbent by vacancy forces. This process is 
much slower than physical adsorption because of the displacement of atoms that 
must occur in the molecules. Chemisorption is an irreversible and endothermic 
process and is usually confined to single layer of adsorbent (Srihari 2007). 

Since most of the natural processes involve surface interaction between 
solid/liquid phase (such as mineralization of magma, formation of fogs and mists, 
cloud condensation, and so forth), and so also the anthropogenic processes (such as 
filtration, separation, coagulation, and so forth), the fractal analysis of concerned 
interaction is fairly justified on account of the self-similarity of these processes as 
well scale independence. In fact, the tagging of surface roughness has been one of 
the important predicaments for both, basic and applied sciences as an effect of the 
fact that many physical and chemical processes take place in porous environments. 
The conventional description of surface irregularity of porous media is based on the 
idea that the chaotic systems come into extensive physical and chemical properties 
of the ordered systems. The deviations in the perfectly ordered systems are very 
small. The researchers have put substantial amount of effort to study the surface 
phenomena analysis in the light of the fractal theory. The main aim of this 
hypothesis is the study of strongly disordered systems or fractal systems, because as 
we go to smaller or larger scales the disorder is more predominant. These systems 
have room for arrangement within structure and occupy inherently more space than 
non-disordered systems. The disorder can be described in terms of a non-integral 
dimension, D, which is a measure of the space-filling ability of the system 
(Zarzycki 1987). The fractal theory has established very successfully, both in its 
application to a wide variety of complex surface geometries and in advancing the 
understanding of how the geometry affects the physical and chemical properties of 
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systems (Pfeifer and Avnir 1983; Pfeifer 1988). In view of the fact that the first 
investigation of the fractal surface properties of solids at molecular scales (Pfeifer 
and Avnir 1983), experimental investigations have covered a wide variety of 
materials with a well-defined fractal surface characteristic as determined by dif- 
ferent techniques (Peleg et al. 1984). Among the properties sensitive to surface 
characteristics, the absorption of gases and vapours can be mentioned. To study the 
effect of the fractal surface on physical adsorption has been grouped into two: one 
has a topological character (Cole et al. 1986) whereas the other, according to Fripiat 
et al. (1986), is based on a molecular approach. Fripiat et al. (1986) and Pfeifer et al. 
(1989) developed a fractal isotherm equation and compared results with the iso- 
therm derived by Cole et al. (1986). They found certain discrepancies particularly at 
high relative pressures. Such a discrepancy has been attributed by Pfeifer et al. 
(1989) to ‘multiple-wall’ effects, neglected in the development of the model. As 
Pfeifer et al. (1989) pointed out; they treat the filling of pores without considering 
that as a pore is being filled, the film grows from two opposite walls and stops 
growing when the two films meet. 

However, conventional approach for estimation of FD does not show the dis- 
tinctiveness of the adsorption process, nor the adsorbents (as established in this 
paper), rather a need for differential adjacency-based fractal analysis would be more 
useful to capture the signature of the regional interactions. 


4 Experimental 


Two distinct agro-based adsorbents, namely, black gram husk (BGH) and green 
gram husk (GGH), were dried, crushed, sieved and desiccated. Phenol, which is a 
toxic chemical, commonly found in many industrial effluents (coke-oven, leather, 
textile, pharmaceutical, etc.), was used as a simulated aqueous adsorbate for the 
study. The study was carried out at 100 ppm phenol concentration for 24 h for 
various adsorption dosages, and different operational conditions (pH, contact time, 
particle size, etc.). Only the cases with 100 % adsorption were used for 
SEM-analysis, with respect to both the adsorbents (i.e., before and after adsorption). 

The samples collected before (BHG-B; GGH-B) and after adsorption (BGH-A; 
GGH-A) were used for SEM imageries. SEM-images consist of pixels where the 
intensity (colour) of each pixel is proportional to the number of back-scattered 
electrons, emitted from the corresponding point at the surface of a material. Some of 
them represent the corresponding points of the micro-particles (foreground pixels), 
others represent background. In the present study, images with 256 grey levels were 
used. For the extraction of the shape of particles only foreground pixels are of 
interest. Two distinct approaches have been adopted to study the FD of the SEM 
imageries. 
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4.1 Evaluation of Total Fractal Dimension of the Imageries 


Using the SEM imageries as population, six systematic samplings were taken (equal 
size) so as to have better representative analysis (Figs. | and 2). Then, using HarFA 
5.5 software (Harmonic and Fractal Image Analysis), the range analysis was carried 
out for each image to evaluate the threshold resolution, and then using this 
threshold resolution as base, the fractal analysis using box-counting method was 
carried out (logarithmic spectrum, mesh size: 0-35, with num steps 30). The FDs 
obtained for respective sub-images were studied for their normality using Shapiro— 
Wilk test as followed in Sect. 4.2 (and suitable comparison tests were adopted). 
The Shapiro—Wilk Test (Shapiro and Wilk 1965) is a test performed to estimate 
the normality of a data set where the test statistics is ‘W’ and null hypothesis may 
not be accepted, for the predefined level of significance, if W is less than the 
corresponding threshold. The test statistics W may defined as follows: 
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Fig. 2. Sub-images of the original SEM imageries of GGH adsorbent 
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sampled from the standard normal distribution and V is the covariance matrix of 
those order statistics. 

Alpha (i.e., the level of significance) for normality has been taken to be 0.05 (i.e. 
the confidence interval being 95 %) and hence when the p-value (the probability for 
the test statistics) is less than 0.05, the null hypothesis has to be rejected (which 
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means the data set is not normal). This test was carried out to estimate normality so 
that the suitable tests for the comparison of the means could be ascertained. In case, 
any one of the data sets used in ‘test for comparison’ is non-normal, Mann—Whitney 
test for two independent samples was carried out, instead of t-test which is restricted 
to only to data sets which are normally distributed. Even while carrying out f-test, the 
quality of variance of the data sets was estimated to use suitable type of the f-test. 


4.2 Evaluation of F-DAS Scores of the Imagery 


The original image was used to estimate the grey levels of the pixels, using 
MATLAB. For each of the original image, these matrices were computed and then 
used to simulate box-counting approach, wherein every 2 x 2 sub-matrices were 
replaced by the averages of the values of the contributing elements of the 
sub-matrices. The same approach was taken again by using 3 x 3 sub-matrices and 
4 x 4 sub-matrices. Now using the corresponding elemental values of these matrices 
thus generated, the negative slope of their logarithms versus the dimensions of 
replaced cells, by least count method so as to evaluate the fractal adjacency scores 
among 2 x 2 and 3 x 3 replaced matrices as well as 3 x 3 and 4 x 4 replaced 
matrices. The difference between these two matrices forms the F-DAS scores. The 
original SEM imageries and the F-DAS scores are presented in Figs. 3 and 4. 
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Fig. 3. F-DAS score contours plots for BGH before and after adsorption 
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Original SEM Image F-DAS Score-Contour Plot 
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Fig. 4 F-DAS score contours plots for GGH before and after adsorption 


4.3 Statistical Analysis of Fractal Dimension 
(FD) and F-DAS Scores 


The FD scores for each of the populations (which include six sampling sub-images) 
were studied for their normality. Accordingly, a comparative test was carried out for 
each of the population types, i.e. BGH before adsorption, BGH after adsorption, 
GGH before adsorption and GGH after adsorption to evaluate the significant dif- 
ferences due to adsorption (Tables | and 2). With regard to F-DAS scores, summary 
Statistics was analysed to evaluate the variation of these scores due to adsorption 
(Table 3). 


4.4 FTIR Analysis of the Adsorbents 


In order to confirm the homogeneity of the adsorbents used with regard to their 
functional groups and compare the changes caused by adsorptions with respect to 
the original adsorbents, Fourier Transform Infrared Spectroscopy (FTIR) was car- 
ried out. In fact, direct information on the presence of surface functional groups can 
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Table 1 Test for normality of fractal dimensions of samples of the adsorbent studied 


BGH-B BGH-A GGH-B GGH-A 
Shapiro—Wilk normality test (D-BW) 
Ww 0.965478 0.916347 0.898424 0.902966116 
p-value 0.860794 0.47945 0.364696 0.391759229 
Alpha 0.05 0.05 0.05 0.05 
Normal Yes Yes Yes Yes 
Shapiro—Wilk normality test (D-BBW) 
Ww 0.560907 0.884109 0.83834 0.88003379 1 
p-value 0.00015 0.288461 0.126294 0.269212334 
Alpha 0.05 0.05 0.05 0.05 
Normal No Yes Yes Yes 
Ww 0.911343 0.919642 0.918197 0.912255497 
p-value 0.445289 0.502789 0.492473 0.45 140039 
Alpha 0.05 0.05 0.05 0.05 
Normal Yes Yes Yes Yes 


Table 2 Test for comparison between the fractal dimensions before and after adsorption for the 


samples 


BGH-B and BGH-A 


Tests for comparison of samples (D-BW) 


GGH-B and GGH-A 


Test t-test (equal variance)-two-tailed 

t-critical 2.228138852 2.228 138852 
p-value 0.001776997 0.715393366 
Significance Yes No 

Tests for comparison of samples (D-BBW) 

Test Mann-Whitney test for two independent samples 

U-critical 5.26002884 5.281446128 
p-value 0.054663936 0.377642391 
Significance No No 

Tests for comparison of samples (D-WBW) 

Test t-test (unequal variance)-two-tailed 

t-critical 2.228 138852 2.364624252 
p-value 0.287974327 0.056797659 
Significance No No 


be obtained from IR studies. IR spectra of adsorbents indicate the possible presence 
of oxygen containing groups, such as carbocyclic, quinone, ether, phenolic, and 
lactone, C = C of aromatic rings, and nitrogen containing groups like pyridine, 
nitrile, cyclic amide. 
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5 Results and Discussions 


5.1 Understanding of the Adsorbents Using F-DAS Scores 


The range of F-DAS for pre-adsorbed BGH is 1.2—1.75 and that of GGH is —1.1 to 
0.2. The post-adsorbed range has shown a distinct increase in both the cases, i.e. 
0.8-1.8 and —0.4 to 1.6, respectively, for BGH and GGH. Although it is interesting 
to note that the FTIR spectra of pre-adsorbed BGH and GGH are very much similar 
(Figs. 5, 6 and 7) with all the functional groups lower in case of BGH than GGH, 
yet the ranges of F-DAS Scores do show the differences, the former being entirely 
positive range and the latter extending over both positive and negative ranges. On 
observation of the plots (% counts vs. F-DAS scores range) of total counts as well 
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Fig. 6 FTIR spectra of BGH before and after adsorption 
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Fig. 7 FTIR spectra of GGH before and after adsorption 
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Fig. 8 Variation in total counts at different F-DAS score ranges for BGH before and after 
adsorption 


as their cumulative values (Figs. 8, 9, 10 and 11), we see a clear differences in the 
plots of BGH and GGH, which is unexplored from FTIR analysis. Although both 
pre-adsorbed BGH and GGH are non-normal, yet the mean, standard error, median, 
standard deviation, kurtosis, skewness, range and Inter-Quartile Range (IQR) are 
higher in case of BGH compared to GGH (Table 3). Besides, the BGH showed a 
single distinct peak whereas GGH showed double peak distribution (primary at 
—0.25 and secondary at —0.05 F-DAS scores). 
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Fig. 9 Variation in total counts at different F-DAS score ranges for GGH before and after 
adsorption 
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Fig. 10 Variation in total cumulative counts at different F-DAS score ranges for BGH before and 
after adsorption 
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5.2 Effect of Adsorption on F-DAS Scores 


As per FTIR spectra between adsorbent before and after adsorption, we find the 
lowering of all the functional groups because of adsorption in case of both BGH 
and GGH, higher being in the latter than the former (Figs. 5, 6 and 7). However, 
F-DAS scores analysis gives a better understanding of this scenario (Figs. 8, 9, 10 
and 11). In both the cases, there is reduction in the peak percentage of F-DAS 
scores, but in case of BGH it is mostly marginal without shifting the mean. 
However, in case of GGH, there is a distinct shift of the mean F-DAS scores from 
negative to positive values effected by the adsorption. The primary and secondary 
peaks of pre-adsorbed GGH showed the reverse trends at 1.45 and 1.25, respec- 
tively. The difference between peak values also showed the reduction by 6.2—2.0 % 
as a result of adsorption. When F-DAS scores were plotted against cumulative 
frequency percentage, we find that the BGH shows a very close distribution pattern 
among pre-adsorbed and post-adsorbed samples, whereas in case of GGH, the 
distinct shifting of distribution towards the right (higher values of F-DAS scores), 
although still the pattern of distribution remains similar. In both the cases, we find 
the very same relationship existing between pre-adsorbed and post-adsorbed 
statistics, namely, increase in median and IQR as well as decrease in mean, standard 
error, standard deviations, kurtosis, skewness and range (Table 3). 


5.3 Scope and Limitations for F-DAS Scores Analysis 


When compared to conventional approach for evaluation of adsorption on adsor- 
bent, F-DAS score can enrich the understanding of the process based on FTIR 
analysis. In fact, the basic method of analysis of FTIR is based on the functional 
group present in the specimen, whereas the F-DAS scores refer to the 
scale-independence persistency of the visual elements of the imagery and hence 
they communicate two different aspects of the samples. The purpose of the work is 
to use the visual analysis to classify the pattern of effect of the adsorbate unto the 
adsorbent, which was more distinct between BGH and GGH than the FTIR anal- 
ysis. In fact, as per FTIR analysis, the functional groups of these adsorbents are 
very much similar, whereas F-DAS scores brought out distinct difference in patterns 
in both the cases in the original samples and corresponding adsorbed samples. 
Thus, this technique can be used for better classifications of various substances 
even with similar chemical compositions and better understanding of the influence 
of stimuli (in this case, adsorption) unto them. 

The F-DAS scores analysis is based on the logarithmic relationship among the 
parameters. And it assumes linearity among all the elements of Y matrices. If the 
relationships for all elements of the Y matrices are normal with regard to the 
corresponding element, the results are most effective. In the present analysis only 
four different box sizes have been utilized to estimate the difference between the 
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slopes. However, if higher box sizes are considered, then the differences will be 
more distinct (i.e. it is more persistent). In case of higher number of box sizes, the 
normality may be assumed (as per central limit theorem) or the values may be 
normalized by suitable transformation. The F-DAS scores analysis carried out in the 
present studies referred to the original grey levels of individual pixels, which may 
be affected by the degrees of exposure during imaging as well as background 
noises. Thus an attempt to use the coefficient of variation scores (i.e. ratio of mean 
and standard deviation) of the image is likely to be more meaningful for F-DAS 
scores studies (as validated by the analysis by authors). The F-DAS scores studies 
may be applied to various samples from original pictures so as to get a better 
understanding of the adsorption phenomena to derive the estimates and standard 
errors at given levels of significance. Furthermore, the segregation method adopted 
in this study for the various pixels enclosed by the selected box size has ‘arithmetic 
mean’, which can be changed to other methods of aggregation such as ‘geometry 
mean’ or ‘centroid’ and so forth, which may give different (even more relevant) 
insights into the domain of application under concern. 


6 Conclusions 


As evident in the preceding paragraphs, we find a distinct advantage of F-DAS 
scores analysis compared to conventional fractal dimensional estimation approach 
(using box-counting method). In case of conventional approach, the FDs of D-BW, 
D-BBW and D—-WBW show different distribution patterns among the samples and 
which are not indicative of the nature of modifications due to stimuli (in our case, 
adsorption). Herein, except D—-BW, the other FDs (i.e. D-BBW and D-WBW) do 
not even show significant variations (at 0.05 level of significance), before and after 
adsorption. On the other hand, F-DAS scores approach provides a better under- 
standing of surface complexities and persistency and the distribution over the range 
of entire spectra. Therefore, although the statistical similarity is being maintained 
by pre- and post-treated (here, adsorbed) samples, they do show distinct variability 
with regard to mean, standard error, median, standard deviation, kurtosis, skewness, 
range and IQR. This variability may be ascribed for specific F-DAS scores sig- 
nature for different adsorbent—adsorbate systems. In fact, the various systems may 
be classified using F-DAS scores approach. The application of these techniques 
may be employed in case of various physico-chemical or biological interactions in 
various domains. 

Here, F-DAS analysis was used as an objective analytic tool to quantitatively 
estimate the behavior of different materials to different responses (in this case, 
adsorptions) using SEM imageries. Traditionally, the SEM imageries have been 
used mostly by visual inspections, which are subjective and depend on the inter- 
preters’ experience. However, using F-DAS scores the researchers will not only be 
able to get a visual understanding of the effect of specific phenomenon/stimuli 
(here, adsorption) but also quantify the effects spatially as well as perceive the 


96 A. Das et al. 


zones of persistency, objectively. Additionally, as mentioned in the preceding 
paragraph, this can also be used to classify various substances, of the similar 
compositions, in terms of their persistency and behavioral complexities, unexplored 
by conventional physico-chemical analyses. 
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The Multi-fractal Scaling Behavior 
of Seismograms Based on the Detrended 
Fluctuation Analysis 


Simanchal Padhy 


Abstract The multi-fractal scaling properties of seismograms are investigated in 
order to quantify the complexity associated with high-frequency seismic signals. The 
third-order MDFA (MDFA3) method is capable of characterising the multi-fractality 
of earthquake records associated with frequency- and scale-dependent correlations 
of small and large fluctuations within seismogram. These correlations are related to 
changes in waveform properties and hence are a measure of the heterogeneities of 
the medium at different scales, sensed by direct and converted phases in a seismo- 
gram with different amplitudes and phases. The non-linear dependence of gener- 
alised Hurst and mass exponent with order g confirms the multi-fractal nature of 
earthquake records. Amongst different types of earthquakes analysed, the 
multi-fractal properties are more pronounced for signals with distinct P-, S- and coda 
waves. The multi-fractal singularity spectrum parameters (maximum, asymmetry 
and width) are used to measure the frequency-dependent complexity of seismo- 
grams. The degree of multi-fractality decreases with increasing frequency, and is 
generally more for the time period windowing dominant seismic phases in the 
seismogram. Significant difference in spectrum width between the original record 
and its randomly shuffled surrogates demonstrates that the multi-fractality in 
earthquake records is predominantly due to long-range correlation of small and large 
fluctuations within seismogram, although its origin due to broad probability distri- 
bution cannot be completely ruled out, based on the values of scaling exponent 
(H, * 0.5) and their weak q-dependency for the surrogates. 
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1 Introduction 


A signal X[n] is scale invariant when X{c +n] = c” - X|n], where H is the Hurst 
exponent, a dimensionless estimator for the self-similarity of a time series, and the 
constant c indicates contraction if c > 1, or dilation if c < 1. The exponent H can 
be estimated from the decay rate of power spectrum of time series. However, the 
spectral analysis suffers from the problem of aliasing that occurs for frequencies 
beyond the Nyquist frequency (half of the sampling frequency). Thus, the standard 
spectral approach is not sufficient to describe the true scaling behavior of time series 
at higher frequencies. In such case, fractal analyses can be used to estimate the 
power-law exponent, H to define the scale-invariant structure of time series data. 

The mono-fractal structures are defined by a single power-law exponent and 
assume that the scale invariance is independent of time and space. However, the 
spatial and temporal variations in scale-invariant structure are characterised by a 
spectrum of power-law exponents (multi-fractal structure) rather than a single 
power-law exponent (mono-fractal structure). 

Many time series data including geophysical signals do not exhibit a simple 
mono-fractal scaling behavior (Hu et al. 2001; Kantelhardt et al. 2001); their scaling 
behavior is more complicated, and different scaling exponents are required for 
different parts of time series (Chen et al. 2002). In such case, multi-fractal analysis 
is needed that requires a multitude of scaling exponents for the full description of 
scaling behavior of such records. The multi-fractal scaling in a time series can be 
ascribed to (1) broad probability density function for the values of the time series or 
(ii) different long-range correlations of the small and large fluctuations in the data. 
Recently, Telesca et al. (2015) showed that the origin of multi-fractality of volcanic 
signals recorded during pre- and eruptive phases in Canary Islands is due to 
long-range correlation related to the change in dynamics during volcanic eruption. 

The scaling properties of both temporal and spatial distribution of earthquakes 
have been studied by using fractal methods (Kagan and Knopoff 1980; Main and 
Burton 1984; Smalley et al. 1987; Rundle 1989; Paladin and Vulpiani 1987, Hirata 
and Imoto 1991; Godano et al. 1996; Telesca et al. 2001, 2004; Dimri 2005; Tang 
et al. 2012; Padhy et al. 2014). 

Very recently a few statistical methods have been performed on seismograms to 
investigate their dynamical characteristics. One of these is the Fisher-Shannon 
method that was applied to identify the tsunamigenic character in seismograms of 
very large earthquakes occurred worldwide (Telesca et al. 2013), to have additional 
information on the seismic hazard (Telesca et al. 2014a) or to detect dynamical 
changes correlated with different stress states of the magmatic setting and the 
plumbing system in volcanoes (Telesca et al. 2014b). 

In this study, we analyse seismic waveforms to explain the complexity of 
seismograms in terms of change in dynamics of the system including medium 
heterogeneities, by using multi-fractal methods. 

The earthquake signal, that characterises a complex dynamical system, often 
exhibits long-range correlations (Sornette 2004) of fluctuations in the data. It is 
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important to quantify such long-range correlations for a better understanding of the 
dynamics of the underlying complex systems. Amongst the various methods used 
to investigate the fractal scaling properties of non-stationary time series, the most 
effective is the detrended fluctuation analysis (DFA). The DFA is a technique able 
to determine the fractal scaling properties and long-range correlations in noisy, 
non-stationary time series (Buldyrev et al. 1995, Bunde et al. 2000; Talkner and 
Weber 2000; Ashkenazy et al. 2001). The DFA method was generalised to 
multi-fractal non-stationary time series data by Kantelhardt et al. (2002), referred to 
as multi-fractal DFA (MDFA). Later on this method has been successfully applied 
to study the multi-fractal scaling behavior of many non-stationary time series data 
(Kantelhardt et al. 2003; Telesca et al. 2004, Lan et al. 2008; Shang et al. 2008; 
Yuan et al. 2009). 

In this study, the MDFA method is used to investigate the multi-fractal scaling 
properties in earthquake records by analysing the multi-fractal spectrum. The 
multi-fractal spectrum is obtained with optimal choice of model parameters suited to 
earthquake data analysed here. We compare the width and shape of multi-fractal 
spectrum obtained for different events to characterise the observed complexity of 
seismogram. We examine the dependence of different parts of seismogram, data length 
(lapse time), frequency band, sampling frequency on the nature of the multi-fractal 
spectrum and its properties, and ultimately interpret the variation in width and shape of 
multi-fractal spectra to infer the nature of medium heterogeneities. Finally, we 
investigate the origin of multi-fractality in earthquake records analysed. 


2 Method of Analysis 


The DFA for a one-dimensional data series can be described as follows. Let X(i) be 
a time series of length N (i = 1, 2, 3, ..., N). The trajectory or profile Y(i) of the time 


series is determined by taking the sum of deviation from the mean value, X 
(Kantelhardt et al. 2002; Telesca et al. 2004) ie. 


i 


Y(i)= 5 (X(k)-—X), i= 1,2,3,...,N (1) 


k=1 


The integrated time series is divided into N, = int (W/s) non-overlapping segments 
of equal size s. Since the length of the series may not be an integral multiple of s, an 
unequal and short part (<s) of the profile may be left at the end of the series. In order 
to consider such part, the same procedure is repeated in the backward direction 
starting from the end. Thus, 2N, segments are obtained altogether. Then the local 
trend for each of the 2N, segments is calculated by a least-squares fit of a polynomial 
function. The variance of each segment v in forward direction is given by: 
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F(sy)p=+ {Ml - Ds + - OY, 2) 


where v may vary from | to N, and y,(z) represents the least-squares fit of the 
segment v for profile Y(i). Similarly, the variance of each segment in backward 
direction is given by: 


where v may vary from N, + 1 to 2N,. 
After detrending the time series, an average is performed over all the segments 
(v = 1-2N,) to obtain the gth-order fluctuation function as: 


1 2M 1/q 
——)-|F(s, | (4) 


where in general, the index variable g may take any real number other than zero. 
Since for g — 0, 1/q blows up, F,(s) cannot be obtained by the normal averaging 
procedure described in Eq. (4), instead a logarithmic averaging procedure is applied 
to get: 


2Ns 
Fo(s) = ot ay do mnlF(s, nih rt (5) 


A similar procedure is performed for different values of time scale length s of the 
segment v. If the time series is governed by long-range power-law correlation, F,(s) 
also follows a power-law for large values of s as: 


F,(s) x sia (6) 


The plot of log, F,(s) versus log, s gives a straight line with slope H,, known as the 
generalised Hurst exponent, for a given value of g. The slope H,, for positive q, 
describes the scaling behavior of segments with large fluctuations, while that for 
negative q, describes the scaling behavior with small fluctuations. From the values 
of H, obtained for each q, one can determine the scaling behavior of the function 
F,(s). A mono-fractal time series is characterised by a unique H, for all values of 
q. In cases, when small and large fluctuations scale differently (multi-fractal time 
series), H, strongly depends on q. For positive values of q, the scaling behavior of 
segments with large fluctuations is characterised by a smaller H,. On the contrary, 
for negative values of q, the scaling behavior of segments with small fluctuations is 
characterised by a larger H, (Kantelhardt et al. 2002). 
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The generalised Hurst exponent H, based on the MDFA method is related to the 
classical or global scaling exponent f, as: 


tq = qH, _ 1, (7) 


A mono-fractal series with long-range correlation is characterised by linearly 
dependent ¢, with single Hurst exponent H. Multi-fractal series, however, is char- 
acterised by non-linearly dependent ft, with multiple Hurst exponent (Ashkenazy 
et al. 2003). 

Using the Legendre transformation, one can introduce the spectrum of local 
dimensions (Holder exponents or singularities), also called as multi-fractal spec- 
trum, D, from the scaling exponent, fy as: 


Dg = ghg = ty; (8) 


where hy = ue is the g-order singularity exponent. The multi-fractal spectrum D, 


expresses the dimension of the subset of the series characterised by singularity 
strength h,, and its width denotes the range of exponents. The spectrum can be 
quantitatively expressed by least-squares fitting it to a quadratic function around its 
maximum /g max (Shimizu et al. 2002) as: 


Dg =a: (hg = hegiel + b 2 (Ng —_ Vaiaus) a Cc, (9) 


where c is an additive constant, c = Di Hewing) = 1. The coefficient b indicates the 
asymmetry of the spectrum. The width of the spectrum, w, obtained by extrapo- 
lating the fitted curve to zero, is defined as: 


w = hg — hg, with Dg(hgi) = Dg(hq2) = 9. 


For a mono-fractal series, H, is independent of q. It thus follows from Eqs. (7) 
and (8) that the width of the spectrum is zero for a mono-fractal series. The three 
parameters (/g. max, w and b) characterising Eq. (9) describe the complexity of the 
signal. If the value of parameter hg max 1s low, then the signal is correlated and the 
underlying process loses fine structure, becoming more regular in appearance 
(Shimizu et al. 2002). The multi-fractal width, w, measures the range of fractal 
exponents in the signal; the wider the range, the higher is the degree of 
multi-fractality. The asymmetry parameter, b, informs about the dominance of low 
or high fractal exponents with respect to the other. A right-skewed spectrum rep- 
resents strongly weighted high fractal exponents, corresponding to fine structures, 
while the left-skewed spectrum represents low fractal exponents with a smoothly 
varying appearance (Telesca et al. 2004). 
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2.1 Estimation of the Multi-fractal Spectrum 
of Time Series (Dq vs. Hy) 


The following procedure is adopted to obtain the multi-fractal spectrum (D, vs. hg) 
of a time series (Vjushin et al. 2001; Kantelhardt et al. 2002). 


(1) Computing the root-mean-square (RMS) variation of time series. 

(2) Finding local detrending of the time series. 

(3) Computing multi-fractal detrending, g-order RMS, F,(s) using Eq. (4). The 
scaling function F,(s) is obtained based on the q-order statistical moments, for 
q ranging between —5 and 5 with an increment of 2. 

(4) Computing q-order Hurst exponent (H,) as the slope of scaling function 
F,(s) for each g, and q-order mass exponent (¢,) using Eq. (7). 


(5) Computing q-order singularity exponent, h, (h = st), q-order singularity 


dimension, D, using Eq. (8), and 
(6) Finally obtaining the multi-fractal spectrum as D, versus h, curve. 


2.2 Selection of Parameters 


The parameters, such as scale s, qg-order and trend order m for obtaining the 
multi-fractal spectrum, are selected based on the following criteria. 


2.2.1 Selection of Scale 


The minimum and maximum segment sizes are selected in order to obtain 
numerically stable RMS values and hence to obtain a stable scale function F;,(s). 
The minimum sample size should be considerably larger than the polynomial order 
m to prevent over-fitting of polynomial trend (Kantelhardt et al. 2002). The max- 
imum segment size (i.e. scale) should be small enough in order to have large 
number of segments in the computation of F’,(s) and hence for a stable estimation of 
the multi-fractal spectrum D, (Scafetta et al. 2003). Here, we used the polynomial 
order for detrending the time series as 3 (m = 3). Based on the above criteria, we 
selected minimum and maximum values of scale as 16 and 1024, respectively, for 
computing the spectrum (D, vs. hg). 


2.2.2 Selection of Order q 


Higher value of g (both positive and negative) may induce larger numerical error in 
the tail of the multi-fractal spectrum. The destabilisation of F, at large positive and 
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negative g-orders also depends on the sample size. For example, time series with 
large sample size have multiple segments of small and large fluctuations, that 
stabilisethe F, computation at large negative and positive q-orders. In addition, the 
stability of the computation of the multi-fractal spectrum depends on the difference 
between the segments of largest and smallest fluctuation. Time series with large 
multi-fractal spectrum width have large differences between the segments with the 
smallest and largest fluctuation and, consequently, the computation of F, at smaller 
negative and positive g-orders is not stable. Here, we used the g-order that weights 
the local variations in time series to vary between —5 and 5 to have optimal 
performance of multi-fractal spectrum. 


2.2.3 Selection of Trend Order m 


A higher value of polynomial order, m, for detrending the time series gives rise to a 
more complex shape of the trend. Multi-fractal spectrum for multiple orders may be 
influenced by non-stationary trends in time series. Linear, quadratic, cubic or 
higher-order polynomials can be used in least-squares fit of time series (MDFA1, 
MDFA2, MDFA3, ...). The potential of trend removal in a series depends on the 
order of polynomial used in the fitting procedure (Kantelhardt et al. 2002). Here we 
selected the value of m in the range 1-3 for minimum segment size of 16, as used 
here for the computation of multi-fractal spectrum (D, vs. hg). 


3 Results 


In Figs. 1, 2, 3 and 4, we show the q-order scaling function F,(s), g-order gener- 
alised Hurst exponent H,, q-order mass exponent f, and the resulting multi-fractal 
spectrum to investigate the multi-fractal scaling nature of noise record and three 
types of earthquake records. Figures la—e show, respectively, the time series, the 
scaling functions F’,(s) (circles) obtained with different segment sizes (scaling, s) for 
different g-orders and the corresponding regression lines computed by MDFA, the 
q-order Hurst exponent, the g-order mass exponent and the multi-fractal spectrum 
for the white noise. In addition, we analysed three different types of seismograms 
based on the MDFA method and compared their multi-fractal spectra. These 
seismograms include (i) an event with clear P-, S- and surface waves with long coda 
waves (Event-1) (Fig. 2), (ii) an event with P-, S- and surface waves, but with more 
diffused and slowly decaying coda waves (Event-2) (Fig. 3) and (iii) an event with 
no clear S- and surface waves (Event-3) (Fig. 4). The complexities associated with 
different types of records may be attributed to varying degrees of medium 
heterogeneities, which are characterised by different degrees of multi-fractal 
structures expressed in terms of H,, tz and width and shape of multi-fractal spec- 
trum. For all types of time series records analysed, the slope of the regression of 
scaling function F,(s) for different q (Figs. 1b-4b) (generalised q-order Hurst 
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Fig. 1 MDFA3 results obtained for the white noise record. a White noise time series, b fluctuation 
functions F;(s) (circles) versus scale s, plotted in log 2—log 2 scales, for different orders g and 
corresponding regression lines (solid lines), ¢ Hg functions, d q-order mass exponent f, versus 
order q and e multi-fractal spectrum (plot of generalised fractal dimensions, D, vs. q-order 
singularity exponents, /t) 


exponent, H,) decreases with q (Figs. l1c—4c), although the exact nature of varia- 
tions depends on the type of signal, explaining the multi-fractal structure of the 
data. The difference between the scaling functions for positive and negative q’s is 
more apparent at small scales (small segment sizes) compared to the large ones 
(Fig. 2b). Such a difference is clearly seen for earthquake signals than noise records 
(see Fig. 1b). The small segments are able to distinguish between the periods with 
large and small fluctuations (i.e. positive and negative q’s, respectively). In contrast, 
the large segments include several local periods with small and large fluctuations 
and therefore average out their differences in magnitude. For white noise, the values 
of H, attain almost constant value of 0.5 with very small variations around 0.5 
(Fig. Ic). The constant H, leads to almost linear mass exponent f, (Fig. Id), that 
further leads to almost constant h, and D,, and finally giving rise to small arcs with 
very small width in the multi-fractal spectrum (Fig. le). The non-constant H, values 
(Fig. 2c) leads to non-linear g-dependency of f, (Fig. 2d) that further leads to 
parabolic (upside-down parabola) multi-fractal spectrum with large width (Fig. 2e). 
Almost a similar scaling behavior (g-dependency of H, and t,) is observed for 
signals of type 2 (Event-2, Fig. 3c—d) and 3 (Event-3, Fig. 4c—d), with a difference 
in the shape and width of their multi-fractal spectra (Figs. 3e—4e). The multi-fractal 
spectrum width (hereafter referred to as MS width) is a measure of the deviation 
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Fig. 2 Same as Fig. 1, but for an earthquake time series (Event-1), recorded with a sampling rate 
of 100 Hz, and characterised by clear P-, S-, surface and coda waves 
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Fig. 3. Same as Fig. 1, but for an earthquake time series (Event-2) recorded with a sampling rate 


of 100 Hz, and characterised by P-, S-, and surface waves, but with more diffused and slowly 
decaying coda waves 
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Fig. 4 Same as Fig. 1, but for an earthquake time series (Event-3) recorded with a sampling rate 
of 100 Hz, and characterised by no clear S- and surface waves 


from average fractal structure for segments with large and small fluctuations. The 
value of MS width is larger for signals (Figs. 2e—4e) than noise records (Fig. le), 
that is expected because of the presence of several seismic phases of varying 
waveform properties in an earthquake record than the noise record. These phases 
are a measure of heterogeneities of the medium of different scales. Moreover, 
high-quality signal with good signal-to-noise ratio, such as Event-1, characterised 
by the presence of clear seismic phases with dominant S-wave, represents large 
anomaly in waveform properties, and shows spectrum with MS width larger than 
other signals analysed (Events-2 and 3). The values of MS width for a normal 
seismogram become larger than those for any seismogram with anomalous seismic 
phases. This property is clearly observed on comparing the MS width for normal 
seismogram (MS width of 1.62 for Event-1) with that of anomalous seismograms 
(MS width of 0.95 for Event-2 and 1.1 for Event-3). 

Next, the shape of the multi-fractal spectrum of the seismogram of Event-1 is 
asymmetrical around the maximum /dmax (Fig. 2e). In particular, it is right-skewed, 
indicating a relative dominance of small fluctuations, while in the case of Event-3 
for anomalous seismograms with strongly attenuated S-waves (Fig. 4e), the 
multi-fractal spectrum exhibits a long left tail indicating that the corresponding time 
series becomes insensitive to the local fluctuations represented by small variation 
and its dynamics is dominated by the large fluctuations. 
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4 Discussion 


4.1 Effect of Sampling Frequency and Window Length 
on MS Width 


The sampling frequency is related to the scale length of heterogeneities of the 
medium sampled by the signal. Here, we attempt to characterise these medium 
heterogeneities from the multi-fractal analysis of seismograms for different sam- 
pling frequencies. The sampling frequency should preferably be selected at least a 
magnitude higher than the dominant frequency of the observed time series so as to 
avoid aliasing-related artefact of the spectrum. A high sampling rate is needed for a 
better understanding of the small-scale medium heterogeneities. Assuming a 
dominant frequency of ~ 10 Hz for local to regional earthquakes analysed here, we 
calculated the MS width for sampling frequencies of 25, 50 and 100 Hz, all greater 
than twice the dominant frequency (Nyquist rate = 20 Hz). It is observed from 
Fig. 5 that the MS width decreases with increase in sampling frequency, consis- 
tently with the nature of variations of scaling function Fq(s) (Figs. 2b—4b) with 
scale. 

Figure 6 shows the variation of MS width with the length of the time series starting 
from the pre-event noise. If the MS width for noise is smaller than that for signal, we 
see that, except Event-1, for the other signals the width increases at beginning from 
minimum corresponding to pre-event noise to maximum corresponding to P-wave 
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Fig. 5 Variation of MS width with sampling rates of 25, 50 and 100 Hz for Event-1. Similar 
variation is observed for other events 
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Fig. 6 Variation of MS width with lapse time (time measured from the beginning of record), used 
for computing multi-fractal spectrum for both noise and signals used in this study 


onset, then decreases corresponding to the arrival of S-wave, and finally reaches the 
lowest value corresponding to later coda arrivals of smaller amplitude at long lapse 
times. Such lapse time dependence of MS width suggests that the change in wave- 
form properties (amplitude and phase) is maximum between pre-event noise and 
P-wave than that encountered between P- and S-wave. However, for Event-1, for 
which S-wave energy dominates the entire seismogram, the MS width reaches 
maximum corresponding to the S-wave arrival group, representing the maximum 
waveform change related to S-wave arrival in a seismogram. Thus, the MS width is a 
measure of the degree of waveform change with lapse time in a seismogram, and its 
variation depends on the nature of seismogram. 

For seismograms of short duration with an unknown trend, application of the 
MDFA method for different detrending polynomial orders will help to distinguish 
spurious spikes or transients from true seismic phases in a seismogram. This aspect 
of automatic seismic phase detection in a seismogram with low signal-to-noise 
ratios, that will accurately estimate the arrival times of different phases is an 
important input to travel-time seismic tomography studies, is planned in future 
studies. 
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Fig. 7 Comparison of multi-fractal spectra for a different parts of seismogram windowed around 
P-, S- and coda waves. b Comparison of MS width for different parts of the seismogram 


4.2 Effect of Different Parts of Seismogram 


We analysed different parts of seismogram windowing different phases, such as 
P-, S- and coda waves (Fig. 7a), to examine the nature of heterogeneities sampled 
by each phase and their effects on MS width and shape (Fig. 7b). We found that the 
change in MS width is large corresponding to the P-wave window in a seismogram 
(Fig. 7b), that corresponds to large change in both amplitude and phase of the signal 
against its background. The coda wave window is found to exhibit minimal width, 
and this is because the change in waveform amplitude and phase for coda waves 
against its background is much less than those for P- and S-waves (Fig. 7a). 


4.3 Effect of Filtering 


We examined the effect of filtering on multi-fractal spectrum by examining the MS 
width obtained from analysing the band-pass filtered seismograms in six frequency 
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Fig. 8 Variation of MS width with central frequency after the seismograms from Event-1 are 
band-pass filtered. The bands used for filtering are 0.5—1, 1-2, 2-4, 4-8, 8-16 and 16-32 Hz with 
central frequencies at 0.75, 1.5, 3, 6, 12 and 24 Hz 


bands with central frequencies at 0.75, 1.5, 3, 6, 12 and 24 Hz. It is observed from 
Fig. 8 that the MS width is maximum at lower frequencies (<6 Hz) and gradually 
decreases with increase in frequency. We observed that the MS width is relatively 
minimum for frequencies around 12 Hz, which is close to the dominant frequency 
of 10 Hz, as assumed in this study. 


4.4 Origin of Multi-fractality 


We investigate the origin of multi-fractality in the analysed seismograms by ran- 
domly shuffling the original time series. During the random shuffling, all correla- 
tions in the time series are destroyed. Hence, if multi-fractality is due to long-range 
correlation, the shuffled series exhibits mono-fractal scaling. If both the original and 
shuffled series exhibit same scaling, then the multi-fractality is due to broad 
probability density function of the time series. Figure 9a—b shows, respectively, the 
multi-fractal spectra and the scaling exponents of the original time series and those 
of ten surrogates. It is clear that the spectra of surrogates, centred around 0.5, are 
narrower than those of the original time series. These results confirm that the 
multi-fractal properties of the earthquake records are dominantly due to long-range 
correlations. However, because the surrogates show weak q-dependence of the 
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scaling exponents Hg, which are not exactly equal to 0.5, the effect of broad 
probability distribution cannot be ruled out (Kantelhardt et al. 2002). In summary, 
the multi-fractal properties of the earthquake records partly result from long-range 
correlations and partly from broad distributions. 


5 Conclusions 


The multi-fractal DFA was applied in scaling analysis of earthquake seismograms 
to characterise their complexity. The study reveals that the degree of multi-fractality 
is higher for seismograms than noise records, due to the co-existence of seismic 
phases of different waveforms. The large spectrum width is found for the part of 
seismogram containing P-waves, representing a sharp change in the waveform 
properties of signal against the pre-event noise. However, for seismograms with 
anomalous seismic phases, the maximum width coincides with time of occurrences 
of those phases. We could resolve the frequency-dependent scaling properties of the 
seismogram, even of finite length, based on different correlations of small and large 
fluctuations within a seismogram. The minimal and stable multi-fractal spectrum 
width is obtained for band-passed signals with central frequencies close to the 
dominant frequency of the signal. By comparing the multi-fractal spectrum for 
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original series to those for randomly shuffled series, we can attribute the 
multi-fractality for the analysed seismograms to long-range correlations, although 
the minimal effect of broad probability distribution on its origin cannot be com- 
pletely ruled out. 
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Fractal Methods in the Investigation 
of the Time Dynamics of Fires: 
An Overview 


Luciano Telesca 


Abstract Fires represent one of the most critical issues in the context of natural 
hazards. Yearly, they affect large areas worldwide causing loss of biodiversity, 
decrease in forests, alteration of landscape, soil degradation, increase in greenhouse, 
etc. Most of these fires have anthropic causes, however there are natural factors, 
above all summer drought, that strongly influence fire ignition and spread. The 
investigation of the time dynamics of fires can be carried out considering the fire 
process per se or focusing on some signal whose variability can be informative of a 
fire occurrence. In the first case, fires are described by a stochastic point process, 
whose events are identified by spatial location, occurrence time and size of burned 
area, or amount of loss. In the second case, time-continuous signals are employed to 
reveal indirectly the occurrence of fires; one of the mostly used signals is the satellite 
normalized difference vegetation index (NDVI) that gives information about the 
“health” of vegetation and, thus, is suited to enhance the status of vegetation after a 
fire stress. In both cases, the concept of fractal can be used to qualitatively and 
quantitatively characterize the time dynamics of fires. Fractals are featured by 
power-law statistics, and, if applied to time series, can be a powerful tool to 
investigate their temporal fluctuations, in terms of correlation structures and memory 
phenomena. In the present review we describe fractal methods applied to fire point 
processes and satellite time-continuous signals that are sensitive to fire occurrences. 


1 Introduction 


To quantitatively characterize fire dynamics, methodologies, which allow to extract 
robust features hidden in their complex time fluctuations, have to be employed. 
Fractality is one of the aspects of such complexity. A fractal is defined as an object 
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whose sample path included within some radius scales with the size of the radius. 
From such definition, clearly fractal processes evidence a scaling behavior, which 
implies that the statistics used to describe them are power laws. In fact, consider a 
statistics f(x), which depends continuously on the scale x, over which the mea- 
surements are taken. Suppose that changing the scale x by a factor a will effectively 
scale the statistics f(x) by another factor g(a), f(ax) = g(a)f(x). The only nontrivial 
solution for this scaling equation is given by f(x) = bg(x), g(x) = x°, for some 
constants b and c (Thurner et al. 1997, and references therein). Therefore, 
power-law statistics and fractals are very closely related concepts. 

A fractal signal can be investigated in terms of its geometrical characterization as 
self-similar curve; but if the fractality of a time series is studied in order to char- 
acterize its temporal fluctuations, we need to perform second-order fractal mea- 
sures, which convey information concerning its correlation properties. 

The well-known method of the power spectral density represents the standard 
manner to investigate the fractality of a time series. The power spectral density 
furnishes information on the distribution of the power of the time series at various 
frequency bands. Since it is calculated by Fourier transforming the series, it allows 
to identify periodic, multi-periodic, or no-periodic behaviors. If the time series is 
characterized by scaling properties the power spectral density is a power-law 
function of the frequency f, Sif”, with the spectral exponent # measuring the 
type and the strength of the time-correlation structures intrinsic in the time series 
fluctuations (Havlin et al. 1999). If 6 = 0 the fluctuations can be considered purely 
random, typical of white noise processes, characterized by completely uncorrelated 
samples. If £ > 0, the fluctuations are called persistent, meaning that positive 
(negative) increments are very likely followed by positive (negative) increments; 
this is typical of systems governed by positive feedback mechanisms. If 6 < 0, the 
fluctuations are called antipersistent, meaning that positive (negative) increments 
are very likely followed by negative (positive) increments; this is typical of systems 
governed by negative feedback mechanisms (Dimri 2000, 2005a, b). 


2 Fire Sequences 


2.1 Representations of Fire Sequences 


A fire sequence can be considered as a realization of a marked temporal stochastic 
point process, which describes events that occur at some random locations in time 
(Cox and Isham 1980) marked by a quantity related to the burned area. Such 
representation was used in modeling several and diverse point processes like 
earthquakes (Telesca et al. 2009a, b; Telesca 2007a), lightning (Telesca et al. 2008), 
starquakes (Telesca 2005), solar flares (Telesca 2007b), and some human and social 
disasters (Telesca and Lovallo 2006, 2007, 2008). 
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A temporal point process may display a fractal behavior if a number of its 
relevant statistics evidence scaling with related scaling coefficients that indicate the 
represented phenomenon contains clusters of points over a relatively large set of 
timescales (Lowen and Teich 1993, 1995; Teich et al. 1996). 

In order to appropriately analyze the temporal properties of a fire sequence, we 
need a proper representation of the series of fire events; depending on this repre- 
sentation, suitable methods will be applied to disclose its fractal behavior. 

A discrete-time process can be derived from the stochastic point process in two 
equivalent ways: (1) using the interevent time series or (2) forming its relative 
counting process. In the first representation a discrete-time series is formed by the 
tule 7; = t;,; — t;, where t; indicates the time of event numbered by the index i. In the 
second representation, the time axis is divided into equally spaced contiguous 
counting windows of duration T to produce a sequence of counts {N,(7T)}, where 
N,(T) represents the number of events falling into the kth window of duration 
T. The duration T of the window is called counting time or timescale. The latter 
approach considers the fire events as the events of interest and assumes that there is 
an objective clock for the timing of the events. The former approach emphasizes the 
interspike intervals and uses the event number as an index of the time. 

Both representations allow us to use different statistical techniques in order to 
investigate the time behavior of a fire sequence, and, in particular, to identify its 
fractality. 

It should be observed that for point processes fractality is a concept used for 
meaning that a process is characterized by time-clustering behavior, where 
time-clustering indicates a time dynamics that is opposite to a homogenous 
Poissonian one. The clustering behavior of a point process leads to a power-law 
(fractal) shape of some of the statistics, used to describe its properties (Thurner et al. 
1997), and that allows estimating the so-called fractal exponent a (Lowen and Teich 
1995; Teich et al. 1996), whose numerical value is an index of the presence of 
clusterization within the process (Lowen and Teich 1993). If the point process is 
Poissonian, the event occurrence times are uncorrelated; for this memoryless pro- 
cess a © 0. On the other side, a # 0 is typical of point processes with self-similar 
behavior; self-similar meaning that parts of the whole can be made to fit to 
the whole in some way by scaling (Mandelbrot 1982). Thus, we can understand that 
the estimation of the fractal exponent a for a fire sequence plays an important role 
in the general characterization of the mechanism underlying the fire phenomenon. 


2.2. Methods 


Depending on the time structure of the representation of a fire sequence (interevent 
times or counting process), different statistical measures can be defined and 
employed to characterize the temporal distribution of a sequence of fire events. In 
this review, we will describe the coefficient of variation (CV), the detrended fluc- 
tuation analysis (DFA), and the multifractal detrended fluctuation analysis 
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(MF-DFA) for fire sequences represented as sets of interevent times; and the Fano 
Factor (FF), the Allan Factor (AF), and the count-based periodogram (PG) for fire 
sequences represented as counting processes. 


2.2.1 Coefficient of Variation (CV) 


A commonly used measure to evaluate the clustering behavior of a point process is 
the CV, defined as 


Cy = si/(t), (1) 


where <z> is the mean of the interevent times and oa, is the standard deviation: a 
Poissonian process (completely random) has a CV = 1, but a clusterized process is 
characterized by a CV > 1. This coefficient does not give information about the 
timescale ranges where the process can be reliably characterized as a clustered 
process. The CV represents a first indicator of the presence of clustering in a point 
process. It only discriminates between clusterized and Poissonian sequences, but it 
does not convey any information about what timescales are involved in the clus- 
tering behavior (Telesca and Lasaponara 2010). This can be considered a limitation 
of such measure, because a complex phenomenon can be deeply known only if the 
different timescales governing its dynamics are well identified and understood. 
Figure | shows the CV of the interevent times of forest fire sequences recorded in 
the 20 regions of Italy from 1997 to 2003; all the fire sequences appear clusterized 
in time, since their CV is larger than 1; however, it is not possible to derive any 
information about the timescale ranges where the clustering behavior is effective. 


2.2.2 Detrended Fluctuation Analysis (DFA) 


The DFA was developed by Peng et al. (1995), and it allows to investigate the 
power-law correlations of nonstationary signals, whose trends should be well dis- 
tinguished from the intrinsic fluctuations of the system in order to find the correct 
scaling behavior of the fluctuations. DFA is a fractal technique for identifying the 
scaling behavior of time series in the presence of possible trends whose origin and 
shape is not very often known, especially if the time series are observational or 
experimental. The method operates on the time series x(7), where i = 1,2, ..., N and 
N is the length of the series. With x,y. we indicate the mean value 


Xave = So xh) (2) 


The signal is first integrated 
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Fig. 1 Coefficient of variation versus region of the interevent time series of forest fires in Italy. 
Values of CV larger than 1 indicate clustering behavior. As shown by the figure, all the sequences 
have a CV > 1. The CV gives only information about the global clustering behavior of the fire 
series, but does not inform about the timescale range where the clusterization takes place (Telesca 
and Lasaponara 2010) 


k 
y(k) = ys [x() — Xave]: (3) 


Next, the integrated time series is divided into boxes of equal length, n. In each 
box a least-squares line is fit to the data, representing the trend in that box. The 
y coordinate of the straight line segments is denoted by y,(k). Next we detrend 
the integrated time series y(k) by subtracting the local trend y,(k) in each box. The 
root-mean-square fluctuation of this integrated and detrended time series is calcu- 
lated by 


F(n) = |= 2 b® — yal). (4) 


k=1 


Repeating this calculation over all box sizes, we obtain a relationship between F 
(n), that represents the average fluctuation as a function of box size, and the box 
size n. If F(n) behaves as a power-law function of n, data present scaling: 
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F(n) xn’. (5) 


Under these conditions the fluctuations can be described by the scaling exponent 
d, representing the slope of the line fitting log[F(m)] to log(n). For a white noise 
process, d = 0.5. If there are only short-range correlations, the initial slope may be 
different from 0.5 but will approach 0.5 for large window sizes. 0.5 < d < 1.0 
indicates the presence of persistent long-range correlations, meaning that a large 
(compared to the average) value is more likely to be followed by large value and 
vice versa. 0 < d < 0.5 indicates the presence of antipersistent long-range corre- 
lations, meaning that a large (compared to the average) value is more likely to be 
followed by small value and vice versa. d = 1 indicates flicker-noise dynamics, 
typical of systems in a self-organized critical state. d = 1.5 characterizes processes 
with Brownian-like dynamics. 

Figure 2 shows the interevent time series (Fig. 2a) and the DFA curves (Fig. 2b) 
for the city fires of Anshan City (China) (Telesca and Song 2011). The city-fire data 
include 6529 fires occurred between 2000 and 2009. By its population, Anshan is 
the 5th largest city in Northeast China, and the 31st largest city in China, covering 
an area of 9,252 km”, and with a population of 3.4 million. Each fire was marked by 
the estimated loss (in Yuan). The interevent times t of the whole sequence show 
regions with low values of t interspersed with those with high values of 7; this 
indicates a certain heterogeneity of the temporal distribution of fire events, leading 
to the presence of certain clusterization of the city fires. The DFA curve, plotted in 
log—log scales, shows two scaling regions with a crossover at the scale nc * 30; at 
scales lower than nc the process is characterized by weak persistence with scaling 
exponent about 0.68, while for scales larger than nc the process is featured by 
strong persistence with scaling exponent about 1.03 that indicates a flicker-noise (1/ 
J) behavior of the point process modeling the city-fire series. Using the relationship 
between the DFA scaling exponent and the spectral exponent d = (1 + a)/2 (Havlin 
et al. 1999), at small scales appa ¥ 0.36 and at large scales G@pra ~ 1.06. An 
estimate of the timescale corresponding to the crossover scale nc can be obtained 
multiplying the crossover for the mean intervent time <z> ~ 771 min, thus giving a 
crossover timescale of about 17 days. The comparison of the DFA results with 
those obtained on a shuffled series (the shuffling procedure destroys the correlations 
and transforms the series into a purely random one) shows that the two curves (the 
original and the shuffled one) do not overlap and the scaling exponent evaluated for 
the shuffled time series is around 0.5, typical of purely random processes. This 
strengthens the robustness and significance of the scaling exponent estimated by 
using the DFA. 

For loss thresholds Ly, assuming the values 100, 500, 1000, 1500, 2000, 2500, 
5000, and 10* Yuan, the DFA curves, plotted in log—log scales, show very similar 
behavior and indicate that the scaling properties of the city-fire process are invariant 
with the size of the loss (Fig. 3). At higher scales the range of variability of the DFA 
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Fig. 2 a Interevent times of the city fires of Anshan City (Northern China). b DFA curves 
(original and shuffled) for the city fires of Anshan City (Telesca and Song 2011) 


exponents is between 0.57 and 0.99. The double scaling behavior observed for very 
low loss thresholds tends to disappear with the increase of the threshold and to be 
substituted by a single scaling region with lower scaling exponent, suggesting that 
the sequence of large city-fires tends to behave as Poissonian processes. 
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2.2.3. Fano Factor (FF) 

The FF is defined as the variance of the number of events in a specified counting 
time or timescale T divided by the mean number of events in that counting time; 


that is 


BE = <N?(T) > — <N;,(T) >? 
() = <N,(T) > 


(6) 


where <...> indicate the average value. In order to evaluate the presence of scaling, 
the timescale T is varied and a relationship FF(T) ~ T is obtained. 

The FF(T) of a fractal point process with 0 < a<1 varies as a function of counting 
time T as: 


FF(T) = 1+ (F) (7) 


To 


The monotonic power-law increase is representative of the presence of fluctu- 
ations on many timescales (Lowen and Teich 1995). The scaling exponent a is the 
so-called fractal exponent. If a > 0 then the represented phenomenon contains 
clusters of points over a relatively large set of timescales (Lowen and Teich 1993, 
1995; Teich et al. 1996). If a * 0, the city-fire occurrence process is Poissonian and 
the occurrence times are uncorrelated. The crossover timescale Tp is the fractal 
onset time, and marks the lower limit for significant scaling behavior in the FF 
(Teich et al. 1996), so that for T « Tp the clustering property becomes negligible at 
these timescales. Thus, 7 is estimated as the timescale over which the FF increases 
as a power-law function of the timescale 7. FF assumes values near unity for 
Poisson processes. 

Figure 3 shows the FF for fire sequences occurred in southern Italy for counting 
times T of duration of 10 min up to 1/10 of the total period. The FF plots clearly 
indicate the fractal behavior of the fire sequences. The FF increases with linear form 
in log—log scales, and this indicates the presence of correlated structures. The early 
flat behavior can be interpreted as the presence of Poissonian dynamics for short 
timescales T up to approximately 50 min. The scaling behavior presents two 
timescale regimes with a crossover at about 24 h in all the plots. This crossover 
seems more evident in 1997 and 2001 FF plots. The fractal exponent is estimated in 
the longer timescale region, obtaining an estimate of opp ranging from 0.69 to 0.91. 
This result indicates the presence of strong time-correlated structures in the 
sequences of fires. 

The fires are marked by a quantity that for forest fires can be the burned area, for 
city fires can be the loss. The variation of the fractal properties of fires depending on 
the intensity of such quantity becomes crucial if we investigate how the time 
characteristics of the fire process change with the level of the intensity of the fires. 
Figure 4 shows the FF for subsequences of fires with burned area A = Ay, varying 
the value of the threshold burned area A,,. Figure 4a shows the FF curves for At, 
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Fig. 3. DFA curves for the city fires of Anshan City depending on the loss threshold (Telesca and 
Song 2011) 


assuming the values 0, 1, 3, 5, 10, 15, 20, 30, and 50 ha. All the FF curves, plotted 
in log—log scales, show clear linear behavior from the timescale T > 10°° min. The 
exponent Gpp(Am) is a function of the threshold Aj,. Figure 4b shows that app 
decreases with the increase of the threshold area A,, from approximately 0.7 to 0.3, 
with a mean value (+ standard deviation) of 0.5 (4 0.1). This behavior of app 
indicates that the time-clustering degree of the events decreases with the increase of 
the threshold burned area. This also indicates that the fire process tends to behave in 
a Poissonian process with the increase of the burned area. 


2.2.4 Allan Factor (AF) 


Another measure, related to the variability of successive counts, useful to detect the 
event clustering in a point process, is the AF (Allan 1966), defined as the variance 
of successive counts for a specified counting time T divided by twice the mean 
number of events in that counting time 


w= 


This measure reduces the effect of possible nonstationarity of the point process, 
because it is defined in terms of the difference of successive counts (Viswanathan 
et al. 1997). As for the FF, varying the timescale T allows producing a relationship 
between AF(7) and T, useful to detect scaling behavior in the sequence. 
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Fig. 4 FF plots for the years 
1997 (a), 1998 (b), 1999 (ce), 
2000 (d), and 2001 (e) of fire 
sequences in southern Italy 
(Lasaponara et al. 2005) 
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Fig. 4 (continued) 2000 


a,,.=0.77+0.01,f 


T (min) 


The AF of a fractal point process varies with the counting time T with a 
power-law form: 


AF(T) = 14 (7) (9) 


T| 


with 0 <a <3 over a large range of counting times T (Lowen and Teich 1995); T; is 
the fractal onset time for the AF and is estimated as the crossover timescale between 
Poissonian and scaling behaviors. It is related to To by (Teich et al. 1997) 


Tit aoa (10) 


As for FF, AF assumes values near unity for Poisson processes. 
The AF represents the most used measure of the clustering of fire point processes 
and was employed in many study-cases worldwide. 
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Figure 5 shows the AF(T) in log-log scales for the 1997-2003 fire sequence in 
an area of central Italy (Telesca and Lasaponara 2006a, b, c). The analysis of this 
AF curve can be considered as paradigmatic because it reveals several features of 
the fire regime. First of all, the non-Poissonianity of the sequence, since the AF 
curve is not flat for all the timescales T. Then, the time dynamics of the fire 
sequence is characterized by the presence of periodicities, indicated by the drops in 
the AF curve (indicated by arrows in Fig. 5); three periodic components are clearly 
detected at 12, 24 h, and about 3 months. The first two frequencies can be put in 
relation to the semidiurnal and diurnal meteorological cycle linked with the state of 
the dead fuel moisture content of open woodland fuel, in response to the normal 
diurnal weather changes, such as the rise and fall of temperature and solar radiation. 
The third frequency is connected with the typical seasonal variability of the veg- 
etation, in particular grass, which has an important influence on the ignition and 
spreading of fires. Two scaling regions (given by the linear parts of the log-log AF 
curve) exist: one ranging between ~ 1.2 days and ~ 13 days, with scaling exponent 
~ 0.8; the other ranging between ~ 13 days and ~ 42 days, with scaling exponent 
~ 1.55. The presence of two different scaling regions indicates that two different 
mechanisms co-exist (due to mixed natural and anthropogenic causes), and the 
different values of the scaling exponent suggest different strengths of the temporal 
fluctuations of the fire sequence. The presence of scaling in the AF statistics evi- 
dences the self-organization of the fire regime. Fires have a direct impact on veg- 
etation, which, in turn, contributes to the future fire activity (Turner 1898). The 
existence of feedback mechanisms involving both fires and the ecological patterns 
(vegetation type, age, physiognomy, etc.) gives rise to correlation structures and 
memory phenomena, identified by the time-scaling behavior. The vegetation pat- 
terns constrain and at the same time are constrained by the processes that generate 
them. A fire occurring in an area, that was never burned before, creates a pattern of 
burned and unburned vegetation patches, which will influence the occurrence of the 
next fire within a continuous feedback dynamics (Malamud et al. 2005; Gill et al. 
2003). 

Figure 6 shows the AF(7) of the 1992-2007 fire sequence in a fire vulnerable 
area of Patagonia (Argentina) (Ghermandi et al. 2008). The AF pattern is even more 
complex showing three time-scaling regions (given by the linear parts of the 
log-log AF curve): the first ranging between ~1 day and ~6 days, with 
scaling exponent ~ 0.5; the second ranging between ~6 days and ~ | month, with 
scaling exponent ~ 1.1; the third ranging between ~1 month and ~4 months, 
with scaling exponent ~ 1.7. The crossover timescales of 6 days and 1 month could 
have probably an anthropogenic nature. The crossover at about 4 months could be 
connected with the mean duration of the fire season (occurring in the austral 
summer). The time dynamics of the Patagonia fire sequence is characterized by the 
presence of periodicities, shown by the drops in the AF curve at about % year 
(clearly connected with the typical seasonal variability, due to cumulative effects of 
weather, climate as well as normal bio-physiological life cycles), and at about 
1.1 years (related to the yearly cycle of vegetation, connected with the meteoro- 
logical and climatic yearly cycles). 
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Fig. 5 Variation of the FF of the fire sequence in Gargano region (southern Italy) with the 
threshold burned area A;,. a Log—log plot of the FF curves varying the threshold burned area Aj; 
all the FF curves show clear linear behavior from the timescale T > 10°° min. b Dependence of the 
scaling exponent app on the threshold burned area A,,; @pp decreases with the increase of the 
threshold area A,, from approximately 0.7 to 0.3, with a mean value (+ standard deviation) of 0.5 
(+ 0.1). (Telesca et al. 2005) 
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Fig. 6 AF of the fire sequence in Tuscany region (central Italy) from 1997 to 2003 (Telesca and 
Lasaponara 2006a, b, c) 


Figure 7 shows the AF(T) for a very long record of forest fires occurred in Ticino 


(Switzerland) from 1969 to 2008 (Telesca et al. 2010). Different time regimes 
co-exist: Poissonian behavior at small and large timescales (the AF is almost 
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Fig. 7 AF of the fire sequence in Patagonia (Argentina) from 1992 to 2007 (Ghermandi et al. 
2008). © 2008 Worldscientific 
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Fig. 8 AF of the fire sequence in Ticino (Switzerland) from 1969 to 2008 (Telesca et al. 2010) 


constant), two scaling regions with scaling exponents about 0.32 between 1 and 
8 days, and about 0.8 between 8 days and 1/2 year. The crossover at about 8 days 
could be related to anthropic activities (Wang et al. 2010). Then daily and annual 
periodicities, indicated by the two drops in the AF curve, are also typical features of 
the sequence of forest fires. 

The correlation between the occurrence of fires and the meteorological condi- 
tions was investigated by using the AF for both the point process of the fires and the 
point process of the relative humidity, the last one obtained after fixing a threshold 
and considering all humidity values below a certain threshold. Wang et al. (2010) 
compared the AF(7) of the relative humidity point process (with a threshold of 
60 %) measured at Miyakonojo (Japan) between 1998 and 1999 (Figs. 8 and 9) with 
that of the forest fire point process recorded in Japan (Fig. 10) (from Song et al. 
(2005)), a fire in Japan has a good probability to occur if the average humidity is 
smaller than 60 %). 

It can be observed that the humidity (Fig. 9) has the same scaling properties 
shown by fires (Fig. 10), with, moreover, the presence of the daily cycle. This result 
validates the existence of a strong link between weather and fires, and it is very 
likely that the daily humidity cycle determines the daily fire periodicity. 


2.2.5 Count-Based Periodogram (PG) 


The count-based PG (Teich et al. 1997), that is the PG of the sequence of the 
counts, is another statistical measure that allows estimating the fractal exponent a. 
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Fig. 9 Relative humidity recorded at Miyakonojo (Japan) from 1998 to 1999 (Wang et al. 2010) 


This measure performs an estimate of the power spectral density, which gives 
information on how the power of the process is concentrated at various frequency 
bands. The calculation of the PG by means of a count-based approach implies the 
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Fig. 10 AF of the forest fire sequence from 1998 to 1999 (Wang et al. 2010) 
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division of the total observation period into a number N of nonoverlapping win- 
dows of length Ty. The sequence of the counts is then Fourier transformed and 
the PG is obtained by the squares of the coefficients of the Fourier representation of 
the series of counts. Generally, in order to use the Fast Fourier Transform (FFT) the 
number N of the windows of duration Ty is a power of 2. For point processes with 
scaling properties the PG decreases as a power-law function of the frequency f over 
a significant range of frequencies, S(Jof“. Of course, for a finite size fractal 
process the power spectral density behaves as an inverse power-law function in a 
limited range of frequencies, approaching an asymptotic value at high frequencies, 
at which the behavior of the process can be considered Poissonian. The numerical 
value of a is an indicator of the presence of clusterization in the process (Thurner 
et al. 1997). If the point process is Poissonian, the occurrence times are uncorre- 
lated; for this memoryless process, the PG is approximately flat for any frequency 
bands and a ¥ 0. On the other side, a # 0 is typical of point processes with scaling 
behavior. This method introduces a bias at higher frequencies, since the fine time 
resolution information is lost as a result of the minimum count-window size. But, 
the estimation of the exponent a principally involves lower frequencies where this 
bias is negligible. 

Figure 11 shows an example of the application of PG measure for the point 
process of Anshan city (China) fires. Each plot shows the AF curves for the original 
sequence (black curve) and for a random point process (blue curve) sharing with the 
original one the same probability density function of the interevent times. The 
random sequence is obtained shuffling the interevent time series of the original 
sequence and then forming the corresponding point process. Each AF curve is 
obtained varying the parameter Ty (and hence, the number of n overlapping win- 
dows N) from 9858 to 615 min. It is evident that if the scaling exponent of shuffled 
sequence does not depend on the length of the window (the scaling exponent is 
approximately around 0.3 in all cases), this does not apply to the original sequence, 
whose scaling exponent decreases from about 0.9 when the window length is the 
largest to about 0.5 when the window length is the smallest. 

This indicates how careful the analysis of the scaling properties of the fire point 
process has to be done when using this measure. 


2.2.6 Multifractal Detrended Fluctuation Analysis 


Multifractals are characterized by high variability on a wide range of temporal 
scales, associated to intermittent fluctuations and long-range power-law correla- 
tions. Many phenomena show multifractality (Currenti et al. 2005; Telesca et al. 
2001). They are generally characterized by sudden bursts of high-frequency fluc- 
tuations, thus evidencing the presence of different scaling behaviors for different 
intensities of fluctuations. The MF-DFA is an extension of the DFA. The only 
modification concerns the fluctuation function, which becomes 
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Fy(s) = 23 [F°(s, vi (11) 


where v = 1,...,Ns, Ns is the number of nonoverlapping windows of duration s and 
the index variable q represents the moment order, whose meaning will be detailed 
below. Varying the timescale s, F,(s) will increase with increasing s. Then ana- 
lyzing log-log plots F,(s) versus s for each value of g, we determine the scaling 
behavior of the fluctuation functions 


F,(s) oc si), (12) 


Monofractal time series are characterized by h(q) independent of q, also called 
generalized Hurst exponent (Kantelhardt et al. 2002). The different scaling of small 
and large fluctuations will yield a significant dependence of h(q) on q. For positive 
q, the segments v with large variance (i.e., large deviation from the corresponding 
fit) will dominate the average F’,(s). Therefore, if q is positive, h(q) describes the 
scaling behavior of the segments with large fluctuations; and generally, large 
fluctuations are characterized by a smaller scaling exponent h(q) for multifractal 
time series. For negative g, the segments v with small variance will dominate the 
average F,,(s). Thus, for negative g values, the scaling exponent h(q) describes the 
scaling behavior of segments with small fluctuations, usually characterized by a 
larger scaling exponent. 

The computation of the multifractal spectrum by means of the Legendre trans- 
form can be performed. The multifractal scaling exponents h(q) defined in Eq. 8 are 
directly related to the scaling exponents t(q) defined by the standard partition 
function multifractal formalism (Kantelhardt et al. 2002) 


t(q) = qh(q) — 1. (13) 


Monofractal series with long-range correlations are characterized by linearly 
dependent g-order exponents 7(q), i.e., the exponents t(qg) of different moments 
q are linearly dependent on q 


t(q) = Hq —-1, (14) 


because h(qg) is independent of g for monofractal series. Long-range correlated 
multifractal signals have a multiple Hurst exponent, i.e., the generalized Hurst 
exponent h(q), 


h(q) = dt/dq 4 const, (15) 


where t(qg) depends nonlinearly on q. 
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The singularity spectrum f(a) is related to t(q) by means of the Legendre 
transform, 


“= — (16) 


f(a) = qa — 1(4); (17) 


where a is the Hélder exponent and f(a) indicates the dimension of the subset of the 
series that is characterized by a. The singularity spectrum quantifies in detail the 
long-range correlation properties of a time series. The multifractal spectrum gives 
information about the relative importance of various fractal exponents present in the 
series. In particular, the width of the spectrum indicates the range of present 
exponents. The width of the spectrum W is a measure of how wide the range of 
fractal exponents found in the signal is; and, thus, it measures the degree of mul- 
tifractality of the series. The wider the range of possible fractal exponents, the 
“richer” is the process in structure. 

In order to understand the type of multifractality underlying the g-dependence of 
h(q), the random shuffle method can be applied. Generally, two different types of 
multifractality in time series can be discriminated: (i) due to a broad probability 
density function and (ii) due to different long-range correlations for small and large 
fluctuations. In the shuffling procedure the values are put into random order, and 
although all correlations are destroyed, the probability density function remains 
unchanged. Hence the shuffled series coming from multifractals of type (i) will 
exhibit simple random behavior with Nshus(g) = 0.5. While those coming from 
multifractals of type (i) will show h(qg) = hghur(q), since the multifractality depends 
on the probability density (Kantelhardt et al. 2002). If both types of multifractality 
characterize the time series, then the shuffled series will show weaker multifractality 
than the original one. 

In cases of fire point processes, the multifractality can be investigated using the 
representation of the interevent times. Figure 12 shows the h(q)-q relationships for 
the original and shuffled sequence of fires recorded in Anshan city (China) (Telesca 
and Song 2011) for q varying between —10 and 10. The original sequence is 
characterized by a larger multifractality than that of the shuffled one. Furthermore, 
Fig. 13, which shows the Legendre spectrum for both series, reveals that the width 
of the multifractal spectrum for the original series is larger than that of the shuffled 
series. Furthermore, since the h(q)-spectrum of the shuffled series is different from 
that of the original one, the multifractality of the city-fire series is given mainly by 
the different long-range correlations for the small and large interevent fluctuations. 
Moreover, the Legendre spectrum of the original sequence appears left-skewed, 
which indicates an asymmetrical weight between the small and the large fluctua- 
tions of the interevent times, in favor of the large fluctuations that dominate the 
multifractality of the fires. 
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Fig. 13 Legendre spectrum f(a)—a for the original and shuffled city-fire sequence. The red line is 
the parabola fitting the spectrum of the original series (Telesca and Song 2011) 
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3. Fires in Satellite Time Series 


Remote sensing techniques allow to investigate properties and dynamics of 
ecosystems and their interannual variability at different temporal and spatial scales. 
Due to their large coverage, high revisitation frequencies and a reliable consistency, 
the estimations of fire-induced variability in ecosystems can be effectively 
approached by using satellite data. The investigation of fires in satellite time series 
has been carried out mainly in analyzing the vegetation characteristics affected by 
fires. The role of vegetation is well known in terms of its ecological impacts. 
Variations in the composition and distribution of vegetation can arise in response to 
natural hazards and anthropic stresses and represent one of the main sources of 
systematic environmental change on local, regional, or global scale (Crowling et al. 
1996). Fire is one of the most critical stresses for vegetation. The effects of fires on 
soil, plants, landscape, and ecosystems depend on many factors, among them fire 
frequency, fire severity, and plant resistance (Huemmrich et al. 1999). 
Characterizing vegetation post-fire behavior is a crucial issue to estimate the fire 
resilience that is the capability of vegetation to recover after fire. Satellite tech- 
nologies can be profitably employed to study the time dynamics of vegetation 
re-growth after fire disturbance at different temporal fluctuations. Moreover, 
fire-induced dynamical processes in vegetation are very complex, since they affect 
the complex soil—surface—atmosphere interactions, due to feedback mechanisms 
involving human activity, ecological patterns, and different subsystems of climate 
(Telesca and Lasaponara 2008). Therefore, the vegetation patterns constrain fires 
and at the same time are constrained by the fire processes that influence them. 
Therefore, the use of fractal methods to investigate the vegetation characteristics in 
fire-affected and fire unaffected sites or to discriminate vegetation patterns before 
and after fire occurrence has represented an important methodological innovation in 
the recent years. 

The remote sensing of vegetation has been traditionally performed by means of 
several vegetation indices obtained through vegetation spectral properties that 
measure the biomass or vegetative vigor. The vegetation indices operate by con- 
trasting intense chlorophyll pigment absorption in the red against the high reflec- 
tance of leaf mesophyll in the near infrared. The simplest form of vegetation index 
is the ratio between two digital values from these two spectral bands. The most 
widely used index is the NDVI (Campbell 1987). The NDVI informs about the 
plant photosynthetic activity and has been found to be related to the green leaf area 
index and the fraction of photosynthetically active radiation absorbed by vegetation. 
Therefore variations in NDVI values reveal variations in vegetation composition 
and dynamics (Myneni et al. 1996). 
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Fig. 14 Location of the forest test sites, fire unaffected Orsomarso (/) and fire-affected Bocchigliero 
(2) NDVI (left) and Corine land-cover map (right) (Telesca and Lasaponara 2006a, b, c) 


Telesca and Lasaponara (2006a, b, c) applied the DFA method on the decadal 
satellite SPOT-VEGETATION pixel time series of NDVI covering two vegetation 
sites in southern Italy, one affected and the other unaffected by fire, in order to 
detect scaling differences in the time fluctuations of the NDVI. The two sites are 
shown in Fig. 14 and an example of the pixel time series is shown in Fig. 15. 

The NDVI shows a periodical behavior that follows the yearly climatological cycle. 
However, before applying the DFA, itis necessary to remove such seasonal fluctuation 
by removing the decadal mean (the decadal mean is calculated for each decade, e.g., 
first decade of January, by averaging over all years in the record), obtaining the NDVIq 
(Fig. 16). The DFA performed on the time variation of NDVIq of the two pixels 
indicates a value of the scaling exponents larger than 0.5 (Fig. 17) that suggests that the 
temporal fluctuations of both time series are positively correlated or persistent. And 
this suggests that the investigated ecosystems are governed by positive feedback 
mechanisms, which tend to destabilize the system under external forces. The feedback 
mechanisms express a positive circular causality that acts as a growth-generating 
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Fig. 15 NDVI decadal SPOT VEGETATION satellite data of a pixel from 1998 to 2003 for 
a Bocchiglieri site and b Orsomarso site (Telesca and Lasaponara 2006a, b, c) 
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Fig. 16 NDVIq decadal SPOT VEGETATION satellite data of a pixel from 1998 to 2003 for 
a Bocchiglieri site and b Orsomarso site (Telesca and Lasaponara 2006a, b, c) 
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Fig. 17 DFA of the pixel time series shown in Fig. 16 (Telesca and Lasaponara 2006a, b, c) 
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- +. 


Fig. 18 Locations of the investigated burned (red crosses) and unburned (green crosses) test sites 
(Telesca et al. 2007) 


phenomenon and therefore drives unstable patterns. Therefore the vegetational pro- 
cesses possess memory of external shocks, which drive the time dynamics of the 
vegetational covers (Telesca and Lasaponara 2006a, b, c). 

Different types of vegetation covers (forest, shrub-land, and mixed) were ana- 
lyzed by Telesca and Lasaponara (2005), who found that independently of the 
vegetation cover, the scaling exponents for fire-affected pixels range around the 
mean value of ~ 1.14, while those for fire-unaffected pixels vary around the mean 
value of ~0.77; and this indicates that fires play an important role in the temporal 
evolution of vegetation. This suggests the inherent character of the fire-related 
vegetation recovery processes, indicating the existence of positive relation between 
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Fig. 19 Example of DFA performed on two pixel time series of two sites: fire-affected (Zagarise) 
and fire unaffected (Monterotondo) (Telesca et al. 2007) 


the amounts of burned and regenerated biomass. Fires drive more unstable patterns 
in vegetational covers, leading to a more efficient fire-induced vegetation recovery 
process. 

Telesca et al. (2007) extended such a study to the whole Italian territory (Fig. 18) 
where a further standardization was performed on the NDVI data, removing not 
only the decadal mean but also dividing by the decadal standard deviation. An 
example of application of DFA on two pixel NDVI time series is shown in Fig. 19. 

Figure 20 shows the scaling exponents of the pixel time series for the sites 
shown in Fig. 18. The DFA exponents range around the mean value of ~ 1.22 for 
fire-affected sites, while those for fire-unaffected sites vary around the mean value 
of ~0.65. The fire-affected sites show significantly larger exponents than those 
calculated for the fire-unaffected sites. From the ecological point of view, the 
behavior of the exponent indicates the resilience of burned ecosystem (the tendency 
of return to the pre-fire initial conditions). In the burned sites the very strong 
perturbation causes the persistence of signal indicated by the exponent larger than 
those obtained for unburned sites, pointing out to the role played by fires in 
changing its vegetation dynamics (Telesca et al. 2007). 
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Fig. 20 DFA exponents of the pixel time series corresponding to the sites shown in Fig. 18 
(Telesca et al. 2007) 


The behavioral trends in the pre- and post-fire dynamics of vegetation were 
analyzed by Telesca and Lasaponara (2006a, b, c) by applying the DFA to five pre- 
and post-fire subseries of the decadal maximum value composite SPOT-VGT 
NDVI data of a forest site of Ventimiglia (northern Italy), mainly covered by 
Quercus ilex and Quercus pubescens, affected by a fire occurred on September 9, 
2001. Figure 21 shows that the scaling exponents for post-fire subseries range 
around the mean value of ~ 1.46, while those for pre-fire pixels vary around the 
mean value of ~ 1.0. 

The fire increased the persistence of vegetation, and this indicates that more 
unstable patterns have been induced in vegetation dynamics, after fire disturbance. 
The value of the DFA exponent, after the fire occurrence, can be used as a quan- 
titative indicator of the vegetation resilience. The vegetation corresponds to 
Mediterranean Quercus forests, dominated by perennial re-sprouting species that 
recuperate the biomass rapidly post-fire; and due to the disruptive perturbation 
caused by fire disturbance, the persistence of signal indicates by the a > 1.4 in the 
burned sites the strong trend of ecosystem to accumulate biomass rapidly (Telesca 
and Lasaponara 2006a, b, c). 

The recovery capacity is slightly reduced after two consecutive fires. Figure 22 
shows the mean scaling exponent calculated for the three NDVI time series of 
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Fig. 21 DFA exponents of 
the pre- and post-fire pixel 
time series in Ventimiglia site 
(Telesca and Lasaponara 
2006a, b, c) 
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Fig. 21 (continued) (d) 
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pixels not affected by any fire, those included within the two fires, and those after 
the second fire in Bolotana (southern Italy) (Telesca and Lasaponara 2008). It is 
clear the discrimination between the three status of vegetation, with the lowest 
persistent degree for the unburned pixels, the highest persistent degree for the pixel 
time series between the two fires. After the second fire the vegetation lowers its 
resilience, due to a more stable response of vegetation to the second fire stress. 
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Fig. 22 Mean and standard deviation of the DFA exponents calculated for the Bolotana pixels 
(Telesca and Lasaponara 2008) 


4 Conclusions 


Investigating fire regime characteristics can support forest fire management for 
several activities including fire prevision and prevention. The knowledge of the 
temporal dynamics of fires is crucial for understanding lightning causes of forest 
fires, improving fire occurrence prediction and fire management planning, inves- 
tigating the role of fires in landscape processes, land-cover changes and degrada- 
tion, improving predictive models of plant succession following wildfires, etc. 

The relevance of applying temporal fractal analyses to fire distributions or 
time-continuous signals sensitive to fires consists in a better understanding of the 
underlying dynamical mechanisms. 

Fire processes show different scaling regimes, suggesting the existence of dif- 
ferent underlying dynamical mechanisms, with different strengths of the power-law 
fluctuations, indicated by the quantitative value of the scaling exponent, which 
suggests a persistent behavior. 

The different mechanisms are very likely mixed natural and anthropogenic 
causes. The periodicity at about | day observed in several cases of fire distributions 
worldwide can be put in relationship with both the diurnal meteorological cycle and 
the diurnal anthropogenic behavior. Pippen (1999) and Plucinski (1999) found a 
diurnal variation in the state of the dead fuel moisture content of open woodland 
fuel, in response to the normal diurnal weather changes, such as the rise and fall of 
temperature and solar radiation. Even the relative humidity sequence shows a daily 
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cycle, thus being one of the causes of the daily cycle in the fire sequence (Wang 
et al. 2010). The annual frequency that characterizes fire sequences is linked with 
the typical seasonal variability, due to cumulative effects of weather, climate as well 
as normal bio-physiological life cycles. The vegetation, in particular the litter and 
the herb layer in the forest stand understory that have an important influence on the 
ignition and spreading of fires, has a typical annual seasonal growth cycle. 

Persistence found in NDVI signal means that the investigated vegetation system 
is governed by positive feedback mechanisms, which tend to destabilize the system 
under external forces. Within this feedback framework, the concept of persistence is 
very useful in order to characterize the stability/instability properties of vegetation 
dynamics. It was observed that fires increase the persistence of vegetation. This 
indicates that more unstable patterns are induced in vegetation dynamics, after fire 
disturbance. And the value of the scaling exponent, after the fire occurrence, can be 
used as a quantitative indicator of the capability of vegetation recovery or resi- 
lience, indicative of a positive relation between the amounts of burned and 
regenerated biomass. The ability of vegetation to recover after fire disturbance is 
crucial in order to avoid or reduce land degradation. Several studies performed on 
post-fire regeneration in the Mediterranean Basin have found that vegetation 
communities generally tend to exhibit a rapid regeneration after fire disturbances. 
Due to the disruptive perturbation caused by fire disturbance, the persistence of 
NDVI signal indicated by a relatively high value of the scaling exponent in the 
burned sites indicates the strong trend of ecosystem in accumulating biomass 
rapidly. 
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